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1 . 0 SUMMARY 


A User's Manual describing the operation of three computer codes (ADD code, 
PTRAK code, and VAPDIF code) is presented herein. This manual is organized for the 
convenience of the user and contains sections describing the general features of 
the computer codes, the input/output formats, run streams, and sample input cases. 
It is presented in loose leaf format so that changes may be made easily 
as additional capabilities are added to the computer programs. 

This User's Manual constitutes Volume II of the final report under NASA 
Contract NAS3-21269. 
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2.0 INTRODUCTION 


This User's Man jal describes the computer codes used to calculate the operating 
characteristics of premixing-prevaporizing fuel/air mixing passages. The calculation 
procedure utilizes three computer codes: the ADD code which calculates the axisymmetric 
or two-dimensional distributions of velocity, pressure, and temperature of the air 
flow; the PTRAK code which calculates the nonequilibrium heat-up, vaporization, and 
trajectories of the liquid fuel droplets in a three-dimensional flow field; and the 
VAPDIF code which calculates the diffusion of fuel vapor or critical autoignition 
species into a moving air stream. A detailed description of the analytical models 
and numerical procedures used in these computer codes is given in Volume I. The 
User's Manual, Volume, II, contains a description of the operation of the computer 
codes . 

The ADD code was originally developed for NASA Lewis Research Center under 
Contract NAS3-15402 (Ref. 1). Important revisions, including the conformal mapping 
coordinate generator, were developed for the U. S. Army Air Mobility Research and 
Development Laboratory under Contract No. DAAJ02-73-C-0037 (Ref. 2). Further develop- 
ment and improvements to the ADD code were funded by United Technologies Research 
Center and Pratt & Whitney Aircraft, Commercial Products Division. Additional 
improvements, including the incorporation of a two-equation model of turbulence into 
the ADD code, was sponsored by NASA Lewis Research Center under Contract NAS3-21853 
(Ref. 3). A complete description of the ADD code is given in Ref. 4 together with 
a new coordinate generator sponsored under contract PEN3-235. 

The original version of the PTRAK code was developed by United Technologies 
Research Center with its own funding. This version was based on a simple equilibrium 
vaporization model for single-component fuels. The nonequilibrium heat-up and 
vaporization model for a distillate fuel was developed under the current contract, 

NAS 3-21269. 

The VAPDIF code was developed exclusively under the current contract although 
it is largely based on the Contractor's experience in developing three-dimensional 
parabolic forward marching computer codes for predicting the behavior of compressible 
f lows . 

This User's Manual is organized into several sections for the convenience of 
the computer code user. Section 3.0 contains a description of the general features 
of each of the computer codes to provide the user with an overview of the types of 
problems which can be solved. Sections 4.0, 5.0 and 6.0 contain descriptions of the 
operation of each of the computer codes. Sections 3.0 through 6.0 should provide 
the user with sufficient information of setup and solve typical problems within 
the scope of the computer codes. Sections 7.0, 8.0, 9.0 present a more detailed 
description of the PTRAK and VAPDIF codes. A description of the ADD code is given 
in Ref. 4 . 

This User's Manual constitutes Volume II of the final report under Contract 

NAS3-21269 . 
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3.0 GENERAL FEATURES 


3.1 General Features of ADD Code 


Program Languag e 

The Annular Diffuser Deck (ADD) code source program Is written in FORTRAN V 
computer language for use on a UNIVAC 1100/81A computer. Some machine specific 
language, such as PARAMETER and INCLUDE F0RTRAN statements, is used. However, these 
statements may be replaced easily by equivalent code for use on other machines. 
Successful conversions of the code to both IBM and CDC computers have been made 
and these versions of the code are available. The ADD code makes use of a UNIVAC 
routine NTRAN which stores and retrieves large data blocks on disc files; however, 
the ADD code is organized so that NTRAN is easily replaced by the equivalent 
F0RTRAN DEFINE FILE. Finally, it should be noted that the ADD code makes use 
of least squares spline fitting and smoothing subroutines provided by IMSL, Inc. 
which are available at all major computer centers. 

Fluid Properties 

The ADD code can treat any compressible fluid with constant thermodynamic 
properties for the gas constant R and the specific heats Cp and Cv, The molecular 
viscosic, , which is temperature dependent, is estimated using Sutherland's law; 
the molecular thermal conductivity is calculated using a constant value for 
Prandtl number. The viscosity of the fluid at standard conditions and Prandtl 
number are input parameters. If these properties are not specified in the input 
data, the ADD code uses the properties of air at standard conditions. 

Types of Flow That May Be Treated 

The ADD code may be used to treat subsonic compressible laminar or turbulent 
swirling flow in axisymmetric ducts or nonswirling flow in two-dimensional ducts. 
The duct may be annular or two-dimensional with both inner and outer walls; or, 
it may be an axisymmetric duct with only an outer wall. Subsonic flows have been 
treated successfully up to the sonic line. The mixing of hot and cold streams 
has also been analyzed using the ADD code. The code is not applicable to flows 
containing regions of separation or reverse flow. 

Duct Geometry Option (I0PT3) 

The flow through axisymmetric or two-dimensional ducts may be calculated 
provided that the principal flow direction is axial; however, the duct cannot 
contain a right-angle turn. The ADD code is not applicable to flows in ducts with 
discontinuities in flow area that produce regions of flow separation. 
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For convenience, provision is made in the code to analyze flown in straight 
annular ducts (I0PT3*1) or in straight wall, annular diffusers (I0PT3-3) using only 
a few input parameters. For ducts of arbitrary shape (I0PT3-2), the coordinates 
(radii) of the inner and outer walls are specified at JLPTS equally-spaced axial 
stations. To assure that the curve representing the duct contour has continuous 
first and second derivatives, a least-squares spline fitting, smoothing and inter- 
polation procedure is Included in the code. This procedure is used whenever the 
number of streamwise stations (JL) is not equal to JLPTS. 

The specification of the duct geometry must include a straight, annular inlet 
section whose length is at least equal to its height. Two-dimensional ducts are 
treated as annular ducts in which the height of the duct is small compared to the 
radius of the duct. Numerical experiments have shown that, if the height of the 
duct is less than 1/100 of the duct radius, the flow is essentially two-dimensional 
to an accuracy of three decimal places. 

Inlet Flow Options (I0PT1) 

Any arbitrary inlet flow conditions may be specified which are consistent with 
the equations of motion and the turbulence model. Two types of input data are 
required: (1) specification of the inviscid free stream and core flow conditions, 
and (2) specification of the laminar or turbulent boundary layer flow parameters. 
With I0PT1=3, 4, 9, the flow is assumed to be turbulent and with I0PT1=7, 8, 
the flow is assumed to be latrinar. With I0PT1*= 3 or 7, the core flow is calculated 
assuming that the stagnation pressure and stagnation temperature are constant 
across the duct. The input Mach number and swirl angle determine the velocities 
and weight flow, and the static pressure is determined by solving the conservation 
equation for radial momentum. When I0PT1* 4 or 8, the inlet core flow is determined 
by specifying KLL data points for fractional distance Y, stagnation pressure P T (Y), 
static pressure P(Y), swirl angle a(Y), and stagnation temperature T^.(Y). For 
I0PT1-9, the core flow is determined by specifying KLL data points for fractional 
distance Y, streamwise velocity Ug(Y) , stagnation pressure P^(Y), swirl velocity 
U^Y), and stagnation temperature T T (Y). Isentropic flow relations and radial 
momentum conservation equations are used to determine the remaining variables. In 
addition, when I0PT1-4, 8 or 9, the corresponding exit flow data must be provided. 
These data are not required by the calculation procedure but are used only by 
plotting routines which can be used to compare calculated and measured exit flow 
profiles. If the exit plane data are not available, the inlet plane data may be 
repeated. 

The boundary layer velocity and temperature profiles are constructed from 
known analytic solutions using the boundary layer displacement thickness (6 ) and 
a power law (1/n) velocity profile. For laminar boundary layers (10PT1«=7, 8) a 
Blasius profile is assumed. For turbulent flows (I0PT1=3, 4, 9), Cole’s boundary 
layer profile is used with the shape parameter determined from 1/n. 
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In many flow situations, it is often more convenient to specify the weight flow 
rather than velocity or Mach number. For these situations, the user may specify 
the weight flow when using I0PT1« 4 or 8. The static pressure profile is automatically 
adjusted to obtain the required weight flow with the other input variables held fixed. 

It should be noted that the initial plane conditions must satisfy the laws of 
motion and be compatible with the turbulence model. Therefore, the ADD code makes 
many checks on the input data to assure satisfactory starting conditions. As an 
example, the initial plane data are checked to determine if the radial momentum 
conservation equation is satisfied. If it is not satisfied, the input static pressure 
profile is replaced by the static pressure profile calculated from the radial momentum 
equation and a DIAGNOSTIC message is printed. The weight flow calculated from the 
initial plane data is checked to see if it is greater than the choked-flow value. 

If it is greater, the calculation stops anu the value of the choked weight flow 
is printed out. Checks are made to assure that the boundary layer profile can be 
matched to the free stream core flow; the necessary adjustments are made automatically 
and the calculation continues. In all cases where adjustments to the input data 
are made and the calculation continues, a DIAGNOSTIC message is printed. When 
no adjustment is possible or when the flow situation is physically impossible, the 
calculation stops and the user is notified with a DIAGNOSTIC message. A list of 
these DIAGNOSTIC messages is given in Section 4.4 

Grid Selection 


The user may determine the calculation grid using input parameters or the grid 
may be determined automatically. In either case the user must specify the number of 
streamwise stations (JL) and the number of streamlines (KL) . Experience has shown 
that a 50 x 50 mesh is suitable for most problems. Default options exist for both 
the distribution of mesh points in the cross flow direction as specified by the 
mesh distortion parameters DDS and the streamwise step size parameter KDS. In 
selecting the mesh distortion parameters DDS, numerical accuracy requires that a 
sufficient number of mesh points exist in the turbulent sublayer. In practice, 
the first mesh point from the wall should be at Y + = 1.0 and at least 20 mesh points 
should be in the boundary layer. These criteria depend on both the flow 
Reynolds number anc wall friction coefficient. Therefore, if DDS is not specified 
in the input data, a value for DDS is calculated using an algorithm which produces 
good results for most cases. The value for the streamwise step size parameter KDS 
depends on the boundary layer thickness and rate of growth of the boundary layer. 

If KDS is not specified, the code selects a value for KDS between each streamwise 
station using an algorithm which produces satisfactorily results for most cases. 


Print Options (I0PT4) 

The frequency and quantity of output are controlled by the print option I0PT4. 
If I0PT4 > 0, the output consists of the mean flow variables including streamwise 
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velocity U 8 , tangential velocity U^, static pressure P, stagnation pressure P-p, 
stagnation temperature T^, and Mach number M at each streamwlse station for JL sta- 
tions; this printout occurs at every I0PT4th station. If 10PT4 _< -1, additional 
information is printed Including the effective turbulent viscosity and thermal con- 
ductivity, the boundary layer solution in universal coordinates U + (Y + ) , and the 
turbulent kinetic energy distribution; this information is printed every I0PT4th 
station. 

Diagnostics 

The ADD code makes numerous checks during the progress of the calculation. 

If the program is able to remedy a detected problem, a DIAGNOSTIC is printed and 
the calculation continues. If a fatal error is detected, the calculation stops 
and a DIAGNOSTIC notifies the user about the nature and location of the error. A 
complete list of DIAGNOSTICS is given in Section 4.4. 

Calculation Option (I0PT9) 

The calculation c the coordinate system may be stored on a data file and 
retrieved for use in subsequent cases. If I0PT9*1, both the coordinates and fhe 
viscous flowfield are calculated. If I0PT9-2, the coordinate system is calculated 
and stored on file NINE and the calculation stops. If I0PT9*3, the coordinates 
stored on file NINE are recalled and the viscous flowfield is calculated. This 
feature is particularly useful when the user wishes to calculate several flows 
using the same duct geometry. 

Data Files 

Since three separate computer codes are used in the analysis, data is passed 
from one computer code to the other through data files. The ADD code generates 
the coordinate system and stores the results on file NINE. The ADD code also 
generates the viscous flowfield solution and stores the results on file EIGHT. 

Both the PTRACK and VAPDIF codes require the data stored on files EIGHT and NINE. 
The PTRAK code calculates the rate of evaporation of fuel droplets for use as the 
source terms for the solution of the diffusion equation by the VAPDIF code. These 
source terms are stored on file SEVENTEEN. The VAPDIF code uses the data stored on 
all three files. It is recommended that these files be registered and catalogued 
files so that the data may be stored permanently over a period of several weeks. 
Proper use of these files allows the user great flexibility in solving problems. 

cart/Stop Options 

A flow calculation may be started at coordinate station J»I0PTI5 and it may 
be terminated at coordinate system station J*I0PT16. If I0PT15 is not specified, 
it is assigned a value I0PT15-1; if I0PT16 is not specified, it is assigned a 
value I0PT16-JL. The calculation of the flowfield may be continued (or restarted) 
at the JM coordinate station by specifying I0PT17-JM. 
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Turbulence Modala (I0PT12) 

The ADD code is provided with four turbulence models described in 
Volume I. For I0PT12- 0, 1, 2 algebraic turbulence models are used based on 
Prandtl's mixing length theory. For I0PT12-3, a two equation model of turbulence 
is used. Option I0PT12-O uses a turbulence model which is well established for 
equilibrium turbulent flowflelds and is therefore recommended for all calculations. 
The other options (I0PT12* 1, 2, 3) are operation but these models have been 
applied to only a few flowfield situations; the use of these models is not recommended 
at the present time. 

B lade Force Options (I0PT2) (10PT5) (I0PT1O) 

Struts, inlet guide vanes, stators, and rotors are modeled in the ADD code 
as a-priori body forces. Three options exist in the code for calculating these 
forces. If measurements of stagnation pressure P T » swirl angle a, and stagnation 
temperature T^ are available, the blade forces can be calculated from blade element 
theory by setting I0PT2«1. If I0PT5-2, the program uses the inltt/exit flow data 
for T0PT1«4. If I0PT5*1 , separate data must be read for the blade force 
calculation. If I0PT2*3, the blade force is calculated from the flow conditions 
and blade geometry using blade element theory and empirical cascade correlations. 

If I0PT2*4, the blade force is calculated using the distributions of exit air swirl 
angle o^CY) and loss coefficient Zg(Y). 

I0PT1O determines whether the blade is stationary (I0PT1O1, stator) or 
rotating (I0PT1O*2, rotor). 

Global Iteration (I0PT14 ) 

The ADD code can treat small regions of separated and reattached flow 
(a separation bubble) using a global iteration procedure. For these cases, KDS 
must be specified. The first iteration is made with I0PT14-O. The second and 
successive iterations are made by repeating the calculation with I0PT14»1. 

As successive passes (iterations) are made, the solution stored on file EIGHT 
is updated. 


3-5 


R82-915362-40 


3.2 General Features of PTRAK Code 


Program Language 

The Particle Tracking (PTRAK) code is written in the F0RTRAN V computer language 
for use on a UNIVAC 1100/81A computer. Some machine specific language, such as 
PARAMETER and INCLUDE FORTRAN statements, is used. However, these statements may 
be replaced easily by equivalent code for use on other machines. The PTRAK code 
makes use of a UNIVAC routine NTRAN which stores and retrieves large data blocks 
on disc files; the PTRAK code is designed so that NTRAN is easily replaced with 
the equivalent F0RTRAN DEFINE FILE statements by modifying only one subroutine. 

Grid Selection 


The PTRAK axial and radial coordinates and computational grid are calculated 
by the ADD code and stored in file NINE. In addition, the PTRAK code uses the 
viscous flowfield solution generated by the ADD code and stored in file EIGHT. 

All information required by the PTRAK code, such as number of mesh points (JL x KL) , 
the starting station (JFIRST) , and the termination station (JLAST) are also stored 
on coordinate file NINE. Information on the third (4>) coordinate must be input 
to the PTRAK code. This information consists of the number of azimuthal grid 
points (LPHI <_ 50), and the azimuthal step size (A4>) . Finally the number of axial 
steps per ADD code streamwise station (KDS) mast be specified. The appropriate 
value of KDS is determined primarily by how rapidly both the droplet trajectories 
and velocities are changing. At the present time, no algorithm tc alter the axial 
step-size automatically exists in the PTRAK code so that KDS must be set to the 
largest value (smallest step size) required for an accurate solution. 

Boundary Conditions (I0PT3) 

Two types of boundary conditions are possible: a wall boundary condition or 
a periodic boundary condition. At a solid wall, a droplet may either strike 
the wall and remain thereon or it may rebound elastically. Droplets remaining on a 
wall may undergo additional evaporation. Droplet-solid wall interaction boundary 
conditions include: elastic rebound with no vaporization (I0PT3*O, 1); elastic 
rebound with vaporization (10PT3* 2, 3). A periodic boundary condition is 
available for swirling flows in annular ducts Thus if I0PT3* 1 or 3, the duct 
is assumed to have four rectangular solid v >xxs. If I0PT3" 0 or 2, the duct is 
assumed to be an annular duct defined by solid walls and two permeable side-walls with 
periodic boundary conditions such that, if a droplet exits through one side-wall, an 
identical droplet will enter through the opposite side-wall; therefore, 
only an annular segment of the flow needs to be analyzed. 
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Location of Fuel Inlectors 

Each of several fuel injectors may be placed arbitrarily in the duct at any 
axial station by specifyiying ' ts coordinates (r, z, $)• The calculation will 
start at the first computed viscous flowfleld station generated by the ADD code. 
Integration of the droplet equations will not start until the axial station is 
reached at which the first fuel Injector is located; it will continue until all 
of the fuel droplets have evaporated. 

Initial Droplet Condit^n s 

The initial droplet conditions for each class of droplets consist of its 
three initial velocity components, diamete , and temperature; up to 1250 droplet 
classes may be input. The PTRAK code uses a forward marching calculation procedure 
so that all droplet classes must have an initial velocity component in the downstream 
direction. Injection angles are limited, therefore, to values less than 85 deg 
relative to the axial direction. 

Compositon of Fuel Droplets 

Droplets of either single component or multicomponent (distillate) fuels may be 
considered by specifying the appropriate thermodynamic and transport properties 
described in Section 4.4 of Volume I. For droplets of distillate fuel, it is 
necessary to provide both a distillation curve and a Cox chart. 

Droplet/Droplet Collisions (I0PT1, 10PT2) 

Models for droplet shattering (controlled by I0PT1) and droplet-droplet 
collisions and coalescence (controlled by (I0PT2) are available in the program. 

Autoigniticr Model (I0PT7) 

Source terms for the rate of production per unit volume of the critical species 
for autoignition can be calculated for use subsequently in the VAPDIF code. 
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3.3 General Feature* of VAPDIF Code 


Program Language 

The Vapor Diffusion (VAPDIF) code le written in the F0RTRAN V computer language 
for use on a UNI VAC l 100/81 A computer. Some machine specific language, such as 
PARAMKTER and INCLUDE FORTRAN statements, is used; however these statements may be 
replaced easily by the equivalent code for use on other machines. The FTRAK code 
makes use of a UNIVAC routine NTRAN which stores and retrieves large data blocks 
on disc files; the PTRAK code is designed so that NTRAN is easily replaced 
with the equivalent F0RTRAN DEFINE FILE statements by modifying only one sub- 
routine. 

Grid Selection 


The VAPDIF code uses the coordinates and computational grid generated by the 
ADD code and stored in file NINE. In addition it uses the viscous flowfleld solution 
generated by the ADD code and stored in file EIGHT together with the vapor and criti- 
cal species source distributions and boundary conditions Renerated by the PTRAK code 
and stored in various files. All information required to locate data at correspond- 
ing griil points is stored in these files. 

Boundary Condition s 

The boundary conditions are specified completely by the ADD code and PTRAK 
code data files. 

Initial Conditions 


It is assumed that the initial concentration of fuel vapor is xero. However, 
the user may specify a constant mass fraction CZERO (card A) which ia distributed 
uniformly over the entire initial plane. 

R est ar t Option 

The calculation may be started, terminated, or restarted at any axial station 
in the duct. These options are controlled by the input parameters 1ADD, l BEGIN, «wd 
1 END on input card 2. The variable IADD is equal to JFIRST from the ADD code 
calculation. IBEG1N is the first calculation station, not including the initial 
plane. for the VAPDIF code. Normally, the calculation is started with I BEGIN equal 
to 1 ADD + 1. I END is the last calculation station in terms of the ADD code 

coordinates. If 1BEG1N N IADD + 1, the code assumes that the solution has been 
calculated and stored in file ELEVEN for station IADD up to station 1BF.G1N-1. Thus 
by specifying I BEGIN and IKN1) the calculation can be started and stopped at any 
point . 
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Print Options (IPRNTX) 

Since the VAPDIF code solves a three-dimensional problem, very large amounts of 
data are necessary to describe completely the solution. For a maximum of a 100 
streamwise stations, 100 normal (radial), and 50 tangential stations, there are 
5 x 10 5 grid points. At each grid point the concentration (C) , the three 
coordinates (n, s, <J>) , the three metric coefficients, the three physical distances, 
and the three cartesian coordinates must be specified. Therefore the total number 
of data points is 6.5 x 10 6 . Print options (card 3) may be used to limit the data 
printout. In all cases, however, a summary table is given which includes the fuel 
vapor flov iate, fuel air ratio, and mass flow weighted average mass fraction of 
fuel as a function of axial distance. 

Autoignition Model 

Subsequent to the calculation of the three-dimensional fuel vapor distribution, 
the Vapor Diffusion code can be used to determine the distribution of critical 
species and thereby estimate whether autoignition of the fuel vapor-air mixture will 
occur . 
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4.0 OPERATION OF THE ADD CODE 


4.1 Runstream for ADD Code 


It is assumed that: 

1. The program will be executed in TPF$; 

2. The executable (absolute) element is MAPADD; 

3. A catalogued file exists for storing the viscous flowfield 
(and is called EIGHT in this example) ; 

4. A catalogued file exists for storing the flowfield geometry 
(and is called NINE in this example). 

Then the following runstream is sufficient to execute the ADD code. 

@ASG,AX EIGHT. ,D/0/TRK/300000 
@ASG, AX NINE. .D/0/TRK/250000 
@USE 8, EIGHT 
@USE 9, NINE 

@ASG,T 10,D/0/TRK/6000 
@ASG,T 11 , D/ O/TRK/ 50000 
@ASG , T 14 ,D/ O/TRK/ 60000 
@ASG,T 22 , D/ O/TRK/ 300000 
(3XQT MAPADD 
(INPUT CARDS) 

@FREE 8 
@FREE 9 
@FREE 10 
@FREE 11 
0FREE 14 
@FREE 22 
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4.2 Input Format for ADD Code 

The input format for the ADD code is described on the input data coding 
forms which follow. These coding forms are organized with one form per input 
data card. Each form contains the names of the variables, the format, and a 
description of the data. The input option card controls the data that must be 
read. Since not all cards are read, the user should make certain that the 
input data agrees with the input options. 

In general the input data is read as follows: 

Card 1 Title Card 

Card 2 Option Card 

Card 3 Mesh Parameter Card 

Card 4 Duct Geometry Card 

+ data as required by I0PT3 

Card 5 Inlet Flow Card 

+ data as required by I0PT1 

Card b Force Data Card (If I0PT2 f 0) 

+ data as required by I0PT2, I0PT5, I0PT1O 


Card 8 


Slot Flow Data Card (option not available) 



Card 2 OPTION CARD FORMAT (3412) 
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I0PT1 (FL0WIN Option) 

* 3 Inlet flow is computed by specifying data on Card 5. (turbulent flow) 

■ 4 Inlet and exit flow profiles are read from 2*KLL data cards following 

Card 5. Input, fractional distance Y, stagnation pressure Pj(Y), 
static pressure P(Y), swirl angle a(Y), and stagnation temperature 
Tf(Y). (turbulent flow) 

= 7 Inlet flow is computed by specifying data on Card 5 (laminar flow) 

* 8 Same as 4 but for laminar flow , 

■ 9 Same as 4 but: Input fractional distance Y, stagnation pressure 

Pj(Y), streamwise velocity U S (Y), swirl velocity U^(Y) and stagna- 
tion temperature T-p(Y) . (turbulent flow) 

I0PT2 (FORCE Option) 

=0 No blade force 

* 1 Calculate blade force from upstreara/downstream flow data; input 

fractional distance Y, stagnation pressure P T (Y) , static pressure 
P(Y), swirl angle a(Y), and stagnation temperature Tj(Y) 

= 2 Not available 

■ 3 Calculate blade force from cascade correlations 

* 4 Calculate blade force from fractional distance Y, exit flow swirl 

angle (^(Y), and loss coefficient Zg(Y) 

I0PT3 (GDUCT Option) Information follows Card 2 

= 1 Calculate a straight, annular duct. 

= 2 Read co-ordinates of duct. 

= 3 Calculate a straight-walled, annular diffuser. 

I0PT4 (PRINT Option) 

Print solution every I0PT4 station. For example, if I0PT4 * 3, every 
third station will be printed. If I0PT4 * -1, the code provides an 
extended printout; this extended printout includes information about 
the boundary layer profiles and the turbulence model. 
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I0PT5 


I0PT6 

I0PT7 

I0PT8 

I0PT9 


I0PT1O 


(STRUT INPUT Option) 

Strut input data (if I0PT2 ■ 1) used to calculate strut forces 

from experimental data measured upstream and downstream of strut. 

= 1 Read in required profiles. 

= 2 The upstream and downstream strut data cards are identical 
to the inlet and exit flow cards and are not read. 

(STRUT Thickness Effects) 

0 Include strut Torces plus thickness effects 

1 Include strut thickness effects only. 

Not Used 

(PLOT Option) 

= 0 No plots requested. 

= 1 Make CALC0MP plots (not available at LeRC) 

(COORDINATE Option) 

= 0 Make an approximate calculation for both streamlines and 

potential lines — do not save flowfield on disk. Used only 
for I0PT3-1. 

* 1 Make exact calculation of streamlines and potential lines — 
store results on logical unit 9 and complete viscous flow 
calculation. 

= 2 Same as 1 but terminate calculation after coordinate calculations 
are completed. 

= 3 Read geometry from logical unit 9 and use in viscous flow 
calculation . 

(R0T0R Option) 

= 0 No rotors or stators. 

= 1 Stators are present. 

= 2 Rotors are present. 
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I0PT11 


I0PT12 


I0PT13 


I0PT14 


I0PT15 


I0PT16 


I0PT17 


(FL0W Option) 

■ 0 Internal flow. 

» i External flow. 

(TURBULENCE Option) 

«* 0 Use two-layer turbulence model. 

* 1 Use two-layer turbulence model with low Reynolds number 

correction. 

■ 2 Use two-layer turbulence model with ntreamline curvature correction. 

= 3 Use two equation turbulence model (applicable to flows in annular 
diffusers only; i.e., diffusers wich both inner and outer walls). 

(SL0T Option) (Not available) 

* 0 No slot cooling. 

■ 1 Slot cooling. 

(GL0BAL Option) 

= 0 Global iterations not used. 

■ 1 Global iterations used - backward differencing for streamwise 

velocity derivatives in vicinity of separation. 

( JFIRST Option) 

Start flow calculation at station I0PT15 — if omitted, I0PT15 = 1. 

( JLAST Option) 

Stop calculation at station I0PT16 — if omitted, I0PT16 * JL. 

(RESTART Option) 

Restart a previously generated case at station I0PT17. 

NOTE: I0PT9 must be equal to 3 and KDS must be the 

same value as used in previous run (see Card 3). 
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Cards 5b INLET/EXIT FLOW DISTRIBUTION F0RMAT(5E1O. 5) 

I0PT1 = 9 2*KLL CARDS FOLLOWING CARD 5. 
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Cards 6c BLADE THICKNESS CARDS F0RMAT(8E1O. 5) 
CARDS FOLLOWING CARDS 6b IF ISHAPE=4 



Card 7 REFERENCE CONDITIONS F0RMAT(2E1O. 5 , 5Fb . 0, 3E10. 5) 
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4.3 Ou ’.put Description for ADD Code 

The output on each page from the ADD code is largely self-explanatory. A 
general description of the output by page is given below. 

Title Page 


This page presents a list of modifications, dates, and report numbers for 
all changes to the ADD code. 

Input Data Page 

This page presents all input data including all options and default input 

values . 

Inlet Flow Pages 

If I0PT1 =4, 8 or 9. the input data is checked for self-consistency. Both 
input and derived results are printed on these pages. 


Performance Page 

Mean and average quantities of the inlet flow which are frequently used to 
measure or normalize duct performance are printed on this page. 

Wall Conditions Page 

The coordinates of the hub and tip wall, mass flow bleed, and wall temperature 
are printed. For adiabatic walls, the wall temperature is not known before the 
calculation and appears as T^ ' 0 on the printout. 

Wall Geometry Page 

This page prints the wall coordinates, wall curvature, and arc length along 
the walls. 

Cap Average Inviscid Flow Pages 

The ADD code calculates the inviscid flow field and prints the solution for 
each streamwise station as determined by the print option I0PT4. 

G ap Average Viscous Flow Pages 

The ADD code calculates the viscous flow field solution and prints the 
solution for each streamwise station as determined by the print option I0PT4. 

The solution appearing on these pages is stored on file EIGHT. 
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Boundary Layer Coordinate Pages 

When I0PT4 < 0, the velocity and shear stress distribution in universal 
coordinates U"^(Y + ) and t + (Y + ) and the effective turbulent viscosity and thermal 
conductivity distributions predicted by the turbulence model are printed. 

Turbulence Properties Pages 

When I0PT4 < 0, the calculated distribution of turbulence kinetic energy, 
Reynolds stress, Reynolds number of turbulence, and turbulent Richardson number 
are printed. 

Mass Flow Average Page 

At the completion of the calculation, a flow summary is given which includes 
mass flow weighted averages of several variables, pressure recovery coefficient, 
and pressure loss coefficient. 

Wall Surface Conditions Page 

This page presents a summary of wall heat transfer conditions including 
wall friction coefficient, wall temperature, integrated wall area, and integrated 
heat transfer through the wall. 

Wall Radiation Summary Page 

This page presents a summary of data required for wall radiation calculations. 
This output is not applicable for cases in which the ADD code is applied to pre- 
mixing fuel preparation passages. 

Boundary Layer Parameter Pages (Hub and Tip) 

These pages summarize the growth of the boundary layer in terms of displacement 
thickness, momentum thickness, and shape factor. 
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4.4 Diagnostics for ADD Code 

Numerous checks are made during the course of the calculation. If a 
minor error c_eurs, a DIAGNOSTIC message is printed and the calculation 
continues. If a fatal error occurs, a DIAGNOSTIC message is printed and 
the calculation stops. A description of these DIAGNOSTICS is given in this 
section. The DIAGNOSTIC message is always of the form: 

♦♦DIAGNOSTIC NO. XX FOR ANNULAR DIFFUSER DECK** 

where xx refers to one of the errors listed. It should be noted that numerical 
values printed with the DIAGNOSTIC message will be in dimensionless form or in 
English units. 

1) I OPT 3 OUTSIDE RANGE OF ALLOWABLE DUCT OPTIONS 

This error is detected in Subroutine ALTMN. The input option must be between 
1 i I0PT3 f 6. 


2) No solution exists in AMFOR 


This error is detected in Subroutine AMFOR. This subroutine solves the Mach 
number function ^ 

N = M (1*-*^-) / ( 1 ♦ VM 2 ) 

for M given N. The function has a maximum at M ■ 1. Hence 

N( 1) *= [ 2 ( 1 ♦ V ) ] 

Solutions do not exit for values of N > N(l). 

3) MASS FLOW EXCEEDS THE MAXIMUM MASS FLOW POSSIBLE 


This error is detected in Subroutine AMINLT which solves the Mach number 
function 


N = M( 1- 


V- 1 

2 


M 2 ) 


2( Y-l) 


for M given N. This function has a maximum for M - 1 given by 

, — li-i — 

N ( 1 ) = (^1) 2(V-1) 

corresponding to choked flow. 


4) ISHaPE AND I0PT2 ARE NOT CONSISTENT 


This error is detected in subroutine CASC. For blade and strut calculations 
use only I0PT2=3 with any ISHAPE, where 


r: 


I 

1 

I 
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3 1 ISHAPE t 6 

Otherwise, the calculation will stop. 

5) FOR BEST RESULTS ADD A STRAIGHT ANNULAR CHANNEL INLET 

This error is detected in Subroutine C00R1. In the construction of the 
duct coordinates, it is assumed that the inlet has no curvature as shown 
in the figure below. This is not a fatal error because small inlet curvatures 
may be tolerated. For best results add a straight annular section to the 
inlet as shown by the dotted lines in the figure. 



Addition of Straight Annular Channel Inlet 


6) PROGRAM ASSUMES INLET FLOW HAS CURVATURE 

This error is detected in Subroutine C00R1. Same as diagnostic 5. 
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7) WALL CURVATURE IS TOO LARGE AT NATION X. 

This error is detected in Subroutine C00R1 usually if the duct has a discontinuous 
change in wall curvature such as shown in the figure below. 



Discontinuous Change in Wall Curvature 


8) Not Used 

9) GREATER THAN 1. PERCENT NORMAL PRESSURE GRADIENT ERROR RECALCULATE STATIC 
PRESSURE 


This error is detected in Subroutine ERPIN. This subroutine integrates 
the radial momentum equilibrium equation. 






-3 5V ? p 5E 
V dr s B An 



and compares (P - F^) to 
If the error given by 


that computed for the input inlet flow (P 


E - 


P i. P 

LI ' H 


(F, 


P H , l 


is greater than 0.01, the input initial static pressure distribution is 
replaced by the above pressure equation and the inlet flow is recalculated. 


10) Not Used 
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11) MASS FLOW REQUIRED EXCEEDS MAXIMUM MASS FLOW POSSIBLE 

Thi 9 error is defected in Subroutine CKINPT. If it is determined that choked flow 
exist!? in the duct, thia diagnostic will be printed; the weight flow muat be 

i educed . 

12) PRESSURE RISE EXCEEDS PERMISSIBLE PRESSURE RISE 

This error is detected in Subroutine CKINPT and Indicated that the deck cannot 
calculate properly the initial flow prof ilea. Check input for errors. 

13) ITERATION OK BACK PRESSURE CALC. FAILS TO CONVERGE 
This error is detected in Subroutine FINV1S. 

In the calculation of strut forces, it has been assumed that the strut exit 
flow is subsonic and unseparated (i.e., U^ > 0) . If these conditions are 
violated, no solution can be obtained. The calculation will stop. 

14) BOUNDARY LAYER TOO THIN FOR MESH SPACING 

This error is detected In Subroutine Fl.lflWIN. The viscous flow calculation 
requires a finite initial boundary layer thickness. In addition, it requires 
enough mesh points to describe the inlet boundary layer velocity profile. The 
deck assumes arbitrarily that at least five mesh points are required. Thus, 
if this diagnostic occurs, increase the number of mesh points, KL, increase 
the mesh distortion parameter, DOS, or Increase the assumed inlet boundary 
layer thickness. If DDS is input equal to zero, the program automatically 
Hets the mesh distortion parameter to the appropriate value for turbulent 
flow. 

15) TOTAL PRESSURE IS LESS THAN STATIC PRESSURE 

This error is detected in Subroutine K1.0W1N. A check is made on the input 
data for 10PT1 ■ 4 to make sure that P^. > P. 

lb) INPUT DATA NOT IN RADIAL EQUILIBRIUM CORRECTIONS APPLIED TO STATIC PRESSURE 

This error is detected in Subroutine KLflWIN. A check is made of the input 
static pressure data for 10PTI • 4. If the static pressure data are not 
in radial equilibrium, it is assumed that the static pressure data are in error 
and that the other inlet data are correct. Then the static pressure profile 
is computed from 
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XV 


AY 

an 


COS^Qf 
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XP an 



with the ID wall static pressure as a boundary condition. 

17) INPUT DDS MUST BE SPECIFIED 

This error is detected in Subroutine FN0RM. At this time there is no 
algorithm to select automatically the mesh distortion parameter DDS for 
laminar flow. 

18) BLADE DATA ERROR IN CKINPT ROUTINE 

This error is detected in Subroutine CKINPT. Blade data have been input 
incorrectly and must be rearranged with Y increasing. 

19) NO UNIQUE SOLUTION FROM MINVRT 

This error is detected in Subroutine MINVRT. If the matrix used to solve 
for the turbulent flow solution is singular, no solution can be obtained. 
This situation may occur due to numerical truncation errors. 

20) LEADING OR TRAILING EDGE INDEX OF STRUT OUT OF RANGE 

This error is detected in Subroutine SLETE. In order to compute blade 
forces, the strut must be located entirely within the duct length. This 
problem may be eliminated by extending the duct as shown in the figure. 
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21) SLOT INPUT NOT IN INCREASING ORDER 

This error is detected in Subroutine SL0TA. 

The slot input data must be arranged in order of increasing axial distance. 
Check input data. The calculation stops if this error is detected. 

22) CHOKED FLOW IN SLOT NO. 

This error is detected in Subroutine SLTFL0. The slot weight flow is 
determined by the ratio of the stagnation pressure of the slot coolant 
fluid to the local wall static pressure. If this pressure ratio is too 
large the flow may be choked at the slot inlet. The calculation will stop. 

23) BOUNDARY LAYER OVERLAP OR TOO LARGE 

This error is detected in Subroutine FL0WIN. For internal flow, the sum 
of the boundary layer thicknesses on the hub and tip walls must be less than 
the duct inlet height. Check input data. 

24) SET TOTAL TEMPERATURE, PRESSURE, ANGLE TO VALUE AT EDGE OF BOUNDARY LAYER 
- CORRECTIONS APPLIED 

This error is detected in Subroutine FL0WIN. For I0PT1 = 4, the calculated 
boundary layer profiles are matched to input inlet flow profiles. 

A good match requires that the stagnation pressure, P-j., be constant in the 
experimentally determined boundary layer region as shown in the figure (dashed 
line) . 
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25) TRUNCATION ERROR CANNOT BE REDUCED BY STEP SIZE 

This error is detected in Subroutine S0LVI. If the step size parameter 
(KDS) is not specified, it is selected automatically by checking the truncation 
error at each step. When an instability occurs, the program attempts to 
reduce the truncation error by reducing the streamwise step size. If the 
truncation error cannot be .educed below a minimum value, the calculation 
stops. 

26) NUMERICAL INSTABILITY 

This error is detected in Subroutine FC0RCT and Subroutine S0LVI and 
is an indication that the prog, am has calculated negative temperatures or 
pressures. The calculator: stops if this problem occurs. 

27) RHOCX ITERATION DID NO. CONVERGE, ERR = 

This error is detected i Subroutine F0RCE. In determining the blade 
force, an iteratior. scheme is used to determine the downstream static 
pressure. If this iteration fails to converge, this diagnostic is printed 
together with the maximum error found in the iteration. The calculation, 
however, is not terminated. 

28) I0PT3 - 2 OPTION NOT IN USE 

This error is detected in Subroutine F0RCE but this option has been deleted 
from the current version of the ADD code. 

29) SOLUTION REQUIRES REVERSE FLOW, INCREASE WFLOW 

This error is detected in Subroutine CKINPT. For flows with radial pressure 
gradients, there is a minimum weight flow below which reverse flow exists. 

This problem can be corrected by increasing the weight flow. The calcula- 
tion will stop. 

JO) LOAD DOWNSTREAM FLOW DATA CARDS 

This error is detected in subroutine CALINV and indicates that the downstream 
flow data cards, required by I0PT1 = 4 or 9, have not been entered. The 
calculation will stop. 

31) SOLUTION FOR BLADE FORCE DOES NOT EXIST 

This error is detected in subroutine F0RCE. The blade force cannot be 
calculated because no inviscid flow solution can be calculated. (Same 
as DIAGN0STIC 29) The calculation will stop. 
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32) GRADIENT OF METRIC COEFFICIENT - 

FOR BETTER RESULTS ADD STRAIGHT CHANNEL INLET 

This error is detected in Subroutine C00R4. It is assumed that the inlet 
duct has no curvature. To avoid problems, add a straight annular section 
to the inlet. The calculation will continue. 

33) INPUT TOO LARGE FOR COLE'S LAW 
SET N < 


This error is detected in subroutine FLOWIN. 

Cole's friction law requires a certain relationship H = H^ (R e ) such that 
there is an upper bound of n<10. For a solution to exist, 

A * <— (1 - — ) >1.573 

U* H12 


Setting 


Then 


1 

H min = 1 - 1.573! 

< lI e 

U* 

n < 2 

Hmin 


The calculation will stop. 

34) WEIGHT FLOW ITERATION MAY NOT CONVERGE IN SUBROUTINE CKINPT CHECK INPUT DATA. 

This error is detected in Subroutine CKINPT. The weight flow iteration that 
determines the static pressure may not converge if the free stream inviscid 
flow is highly distorted. An input flow which is more uniform in stagnation 
pressure is required. The calculation will stop. 

35) WFLI AND IjOPTll OPTIONS INCOMPATIBLE 

This error is detected in subroutine ALTMN. The weight flow cannot be specified 
for ex .ernal flow. The calculation will stop. 

36) I0PT1 = 1 OR I0PT1 = 2 OPTIONS NOT USED 

This error is detected in Subroutine ALTMN. The options I0PT1 = 1 and I0PT1 = 2 
have been deleted from the code • 
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37) CONFLICT OF OPTIONS, I0PT14<O IMPLIES SEPARATION AND GLOBAL ITERATIONS. 
AUTOMATIC STEP SIZE ALGORITHM CANNOT BE USED. 

This error is detected in Subroutine S0LVI. When performing a global 
iteration (IOPT14>0), the same number of streamwise steps must be used 
for each iteration. Hence the automatic step size algorithm for the stream- 
wise direction must not be used. The calculation will stop. 
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4.5 Debug Options for ADD Code 

When set equal to unity, these options allow intermediate results calculated 
by the subroutine indicated to be printed as an aid in debugging a troublesome 
case. Note that these outputs are not converted to metric units and reference 
must be made to the source code for interpretation of printout. 


OPTION 

SUBROUTINE 

OBJECTIVE OF SUBROUTINE 

IDBGI 

TURB 

Calculates two-layer turbulence model 

IDBG2 

FC0RCT 

Calculates shear stresses and heat fluxes at 
each station 

IDBG3 

FL0WIN 

Generates initial flow profiles 

IDBG4 

SLTFL0 

Calculates slot inlet flows 

IDBG5 

S0LVI 

Calculates viscous flow solutions 

IDBG6 

C00R 

Generates required geometric parameters 

IDBG7 

F0RCE 

Calculates forces generated by struts and 
blades 

IDBG8 

MINVRT 

Inverts a matrix 

IDBG9 

SM00TH 

Smooths duct contour read via I0PT3 = 2 

IDBG10 

GDUCT 

Calculates duct geometry 

IDBGI 1 

SLTFL0 

Obtains additional information from SLTFL0- 
see IDBG4 

IDBG12 

S0LVI 

Obtains additional information from S0LVI - 
see IDBG5 

IDBG13 

CKINPT 

Checks inlet flow input for errors 

IDBG14 

S0LVI 

Debugs the algorithm that automatically 
computes the maximum step size in the stream 
wise direction while assuring computational 
stability . 

I DBG 15 

Not used. 


I DBG 16 

Not used. 


IDBG17 

Not used. 
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4.6 Sample Input for ADD Code 

Two sample inputs to the ADD code are presented on the following pages. 

These cases correspond to the two design studies described in Vol. I Section 7.0. 

The first sample is the input for the Swirl Tube Premixing Passage Case and the 
second sample is the input for the Series Staged Premixing Passage Case. 

Swirl Tube Premixing Passage 

The option card (line 2) indicates that the inlet flow conditions are to be 
calculated from input flow profile data (I0PT1=9). The duct geometry is to be 
determined using the straight wall annular diffuser option (I0PT3=3). Only the 
solution at every 5th station is to be printed (I0PT4=5) and the coordinates are 
to be read from a previously calculated data file (I0PT9=3). The mesh parameter 
card (line 3) indicates that the default mesh distortion parameter (DDS = 0.0) and 
step size (KDS = 0) algorithms are to be used. Two sets of 26 input data cards each 
describing the inlet and exit flow profiles, respectively, will be read. The duct 
geometry card (line 4) indicates that the overall duct length is 11.1 cm and that 
the inlet radius is 2.8956 cm. The inlet flow card (line 5) shows that the initial 
boundary layer displacement thickness is 0.029 cm and that a 1/7 power law profile is 
assumed for each wall. Lines 6 through 31 are the KLL=26 inlet flow profile data 
cards and lines 32 through 57 are the KLL-26 exit flow profile data cards which are 
identical to the inlet flow cards. The last card indicates that default values are 
used for the remaining input parameters exce- c for the Prandtl number, heat capacities, 
and molecular viscosity. 

Series Staged Premixing Passage 

The option card (line 2) indicates that the inlet flow is to be calculated 
assuming a constant stagnation pressure and stagnation temperature in the core flow 
(I0PT1=3). The duct geometry is to be read from input data cards (I0PT3=2). On 
the mesh parameter card (line 3), the default mesh distortion parameter (DDS=0) 
has been selected but the streamwise step size parameter has been input KDS=2. 

From line 3, it is noted that the duct coordinates at JLPTS=50 equally spaced axial 
stations are to be read and that the least squares spline smoothing routine will be 
used (JLPTS i JL) . The length of the duct is 22.72 cm (line 4). Lines 5 through 
11 contain 50 data points for the tip radii and lines 12 through 18 contain 50 
data points for the hub radii. The inlet Mach number is 0.102 (line 19), the 
stagnation pressure is 11.06 atm, and the stagnation temperature is 745K. 
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5.0 OPERATION OF PTRAK CODE 
5.1 Runstreams for PTRAK Code 

Autoignition Model Not Used 
It is assumed that: 

1. The program will be executed in TPF$; 

2. The executable (absolute) element is PTRACK7 ; 

3. The catalogued files, EIGHT and NINE, exist and have been written into 
by the ADD code (see ADD code input format) ; 

4. A catalogued file exists for storing the fuel vaporization terms for use 
with the VAPDIF code (and is called SEVENTEEN in this example) . 

Then the following runstream is sufficient to execute the PTRAK code. 

@ASG,AX EIGHT., D/O/TRK/ 300000 

@ASG, AX NINE., D/0/TRK/250000 

@USE 8, EIGHT 

@USE 9, NINE 

@ASG,T 12.D/700000/TRK 

@ASG , T 14.D/60000/TPJC 

@ASG, AX SEVENTEEN., D/O/TRK/ 300000 

@USE 17, SEVENTEEN 

(?ASG,T 18 

@ASG,T 19 

(3XQT PTRACK7 

(INPUT CARDS) 

@FREE 8 
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3FREE 9 
(?FREE 12 
@FREE 14 
@FREE 17 
@FREE 18 
@FREE 19 

An toigni t ion Model User* 


It Is assumed that: 

1. The program will be executed in TPF$; 

2. The executable (absolute) element is PTRACK7 ; 

3. 'Hie catalogued files, EIGHT and NINE, exist and have been written into 
by the ADD code (see ADD code input format) ; 

4. A catalogued file exists for storing the fuel vaporization terms for use 
with the VAPDIF code (and is called SEVENTEEN in this example) . 

5. A catalogued file exists for storing the critical species source terms 

for use with the VAPDIF code (and is called TWENTYSEVEN in this example); 

6. A catalogued file exists for storing the temperature depression source 

terms for use with the VAPDIF code (and is called TWQJTYEIGHT in 

this example). 

Then the following runstream is sufficient to execute the PTRAK code when using the 
autoignition model. 


@ASG, AX EIGHT., D/0/TRK/300000 
@ASG, AX NINE., D/0/TRK/250000 
@USE 8, EIGHT 
@USE 9, NINE 

@ASG, ; 12 ,D/700000/TRK 
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0ASG.T 14.D/60000/TRK 

9ASG.AX SEVENTEEN. , D/O/TRK/ 300000 

0USE 17, SEVENTEEN 

8ASG.T 18 

8ASG.T 19 

0ASG.AX TWENTYSEVEN. .D/0/TRK/300000 

@USE 27 .TWENTYSEVEN 

0ASG.AX TWENTYEIGHT. .D/O/TRK/ 300000 

@USE 28, TWENTYEIGHT 

8XQT PTRACK7 

(INPUT CARDS) 

8FREE 8 
0FREE 9 
8FREE 12 
8 FREE 14 
8FREE 17 
8 FREE 18 
8 FREE 19 
0FREE 27 
0FREE 28 
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5.2 Input Format for PTRAK Cod* 

The input to the PTRAK code is described on the input deta coding forms which 
follow. These coding forms are arranged with one form per inout data card. Each 
form contains the namos of the input variables, the format, and a description of 
the data. In general, the input data is read as follows: 


Card 1 
Card 2 
Card 3 
Card 4 
Card 5 
Card 6 
Card 7 
Card 8 
Card 9 
Card 10 
Card 11 
Card 12 
Card 13 

Card 14 

Card 15 


Title Card 
Option Card 

Fuel Class Description Card 

Injector Description (IL0C) Cards 

Initial Velocity (IL0C) Cards 

Fuel Flow Rate Card 

Fuel Thermodynamic Constants Card 

Fuel Thermodynamic Functions (6) Cards 

Air Thermodynamic Constants Card 

Air Thermodynamic Functions (3) Cards 

Mesh Description Card 

Distillation Curve (IDSTL+1) Cards 

Collision Data Card 

Cox Chart Data Cards iIC0Xl+l) Cards 

Autoignition Model Constants (3) Cards 


«... 
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Cox Chart for the Simulation of a Pure Substance 


The Cox Chart for simulating a pure substance consists of (a) a set of 
arbitrary "carbon" numbers, C0X(1); (b) a set of arbitrary temperatures for a vapor 
pressure of P1CH, C0X(2); and (c) a set of temperatures all of which are equal to 
the temperature that produces a vapor pressure of P2CH, C0X(3). At any point in 
the calculation, the program will obtain the distillation temperature for the 
instantaneous value of percentage of fuel evaporated (this temperature for a pure 
substance is the normal boiling point), calculate an obviously meaningless "carbon" 
number, and then find the temperature that produces a vapor pressure equal to 
P2CH. These parameters are sufficient to calculate the vapor pressure at the 
instantaneous value of droplet temperature (see Volume I, Section 4.5) • A Cox chart 
for the normal paraffin series is given in Table 5.1. 
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TABLE 5.1 


COX CHART FOR NORMAL PARAFFIN SERIES 


P1CH - 1 atm 
P2CH - 6.8 atm 


Carbon No. 


T(gPlCH (deg K) T@P2CH (deg K) 


2 

175 

3 

230 

4 

270 

5 

310 

6 

345 

7 

370 

8 

405 

9 

430 

10 

455 

11 

475 

12 

490 

13 

510 

14 

530 

15 

545 

16 

560 

17 

580 

18 

588 

19 

605 

20 

620 

23 

640 


225 

285 

335 

380 

420 

455 

485 

515 

540 

570 

590 

605 

625 

640 

660 

675 

690 

700 

715 

740 
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5.3 Output Description for PTRAK Code 

The output on each page from the PTRAK code Is largely self-explanatory. A 
general description of the output by page Is given below. 

Title Page 

This page presents a list of modifications, dates, and references to all changes 
to the PTRAK code. 

Option Page 

This page lists all options used in the PTRAK calculation and the number of 
classes in each category. 

Fuel Injector Page 

This page lists the location and initial conditions for all fuel injectors. The 
total fuel flow rate is also printed. 

Fuel Thermodynamic Properties Page 

This page presents the thermodynamic data and transport properties for both the 
fuel liquid and vapor phases. 

Air Thermodynamic Properties Page 

This page presents both the thermodynamic data and transport properties for the 
air and the mesh parameters. 
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Multicomponent Fuel Properties Page 

This page lists the data for the distillation curve and Cox chart. 

Autoignition Model Constants Page 

This page lists all constants used in the autoiRnition model. 


Initial Conditions by Class Page 

This page describes the initial conditions, including number density, for all 
droplet classes. 

Duct Geometry Page 

This page presents the coordinates of the duct contour. 

Solution Pages 

These pages present all of the dependent variables by class, and the overall 
Sauter mean diameter, number of droplets, and fuel flow rate at each axial station. 

S ummary Page 

The summary page presents the global properties of the fuel spray and includes 
the Sauter mean diameter, liquid fuel flow rate, percentage of fuel evaporated, and 
fuel-air ratio as a function of axial distance. 
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5.4 Diagnostics for PTRAK Code 

Numerous checks are made during the course of the calculation. If a minor 
error occurs, a DIAGNOSTIC message is printed and the calculation continues. If a 
fatal error occurs, a DIAGNOSTIC is printed and the calculation is stopped. A 
description of these errors is given in this section. The DIAGNOSTIC messege is 
always of the form 


**DI AGNOSTIC NO. XX FOR PTRAK CODE** 
where xx refers to one of the errors listed below. 

1. FAILURE TO INTERPOLATE IN FINTP 

This error occurs when the PTRAK code cannot find the location where a particle 
track crosses a grid point. The axial step-size should be reduced by increasing 
KDS. This error causes the calculation to terminate. 
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5.5 Debug Options for PTRAK Code 

When set equal to unity, these options allow intermediate results calculated 
by the subroutine indicated to be printed as an aid in debugging a troublesome 
case. These outputs are not converted to metric units and reference must be made 
to the source code for interpretation of the printout. 


OPTION 

SUBROUTINE 

IDBG1 

PTRAK 

IDBG2 

B0UNCE 


IDBG3 


C0LLSM , COLLDB 
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5.6 Sample Input for PTRAK Code 

This sample of input to the PTRAK code is based on the Series Staged Premixing 
Passage case described in Volume I, Section 7. The option card (line 2) indicates 
that the duct is an annular passage with periodic boundary conditions (lOPT3 a 4). 
Thus the inner and outer boundaries are to be treated as solid walls and the two 
remaining boundaries represent the lateral extent of the annular segment being con- 
sidered. While the air flow in this case is not swirling, droplets injected 
normal to the flow direction but in the angular direction may exit from this seg- 
ment to an adjacent segment. If droplets do exit through one of the "sidewalls", 
then a source of droplets identical in all respects to this sink of droplets must 
enter this segment through the opposite sidewall. Only the behavior of the spray 
in a l/60th segment of the annular duct is to be calculated (LPHI - 31 and DPHI * 
0.2 from line 24). From line 5, the fuel class description card, there are four 
injectors (IL0C = 4) and only one class is to be formed at each injection location 
using the binominal distribution function for each property. Thus, the total num- 
ber of droplet classes is four. These four injectors, described on lines 4 through 
7, are located 7.48 cm from the inlet. Lines 8 through 11 describe the initial 
velocity components and droplet size for each injector. The fuel flow rate is 
.01146 kg/sec (line 12). The thermodynamic and transport data for the fuel and air 
are presented on lines 13 through 23. The distillation curve appears on lines 25 
through 31 and the Cox chart is listed on lines 33 through 53. All droplets which 
strike a solid wall will rebound elastically (line 32, Cl = 1.0). 
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6.0 OPERATION OF VAPDIF CODE 

6.1 Runstreams for VAPDIF Code 
Calculation of Fuel Vapor Mass Fraction Distribution 

It is assumed that: 

1. The program will be executed in TPF$; 

2. The executable (absolute) element is DIFFUSEMAP; 

3. The catalogued files, EIGHT and NINE, exist and have been written into 
by the ADD code; 

4. The catalogued file, SEVF TTEF.'J, exists and has been written into by the 
PTRAK code; 

5. The catalogued file, FUELS^L,, exists and will be written into by the VAPDIF 
code. It is used to permit the VAPDIF code to restart a case subsequent 

to the initial run. If no restart capability is desired, this file may be 
a temporary file. The file FUELS0L contains the solution of the diffusion 
equation at each mesh point for the fuel vapor mass fraction. 

Then the following runstream is sufficient to execute the VAPDIF code. 

@ASG, AX EIGHT, D/O/TRK/ 300000 

MSG , AX NKIE,D/0/TRK/250000 

@USE 8, EIGHT 

@USE 9, NINE 

@ASG, T 10 .D/O/TRK/ 2 2 50000 
@ASG, AX FUELS^L, D/O/TRK/ 2 50000 
0USE 11 , FUELS0L 
@ASG,T 12 /D/O/TRK/ 7 50000 
OASG.T 13, D/O/TRK/ 2 2 50000 
'JASG, AX SEVENTEEN. , D/O/TRK/ 300000 
(3USE 17, SEVENTEEN 
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@XQT DIFFUSEMAP 
(INPUT CARDS) 

@FREE 8 
@FREE 9 
0FREE 10 
(3FREE 11 
@FREE 12 
@FREE 13 
@V REE 17 

Calculation of Critical Species Concentration Distribution 

The mass fraction of critical species (ethene when using autoignition Model I) 
or of the r elat ive mass fraction of the unknown critical species (using autoigni- 
tion. Model II) can be calculated using the following runstream once the fuel vapor 
mass fraction distribution is known. 

It is assumed that: 

1. The program will be executed in TPF$; 

2. The executable (absolute) element is DIFFUSEMAP; 

3. The catalogued files, EIGHT and NINE, exist and have been written into 
by the ADD cede; 

4. The catalogued files, SEVENTEEN, TWENTY SEVEN and TWENTYEIGHT have been 
written into by the PTRAK code; 

5. The catalogued file, ELEVEN, exists and will be written into by the VAPDIF 
code. It is used to permit the VAPDIF code to n.start a case subsequent 
to the initial run. If no restart capability is desired, this file may be 
a temporary file. The file ELEVEN contains the solution to the diffusion 
equation at each mesh point for the concentration of critical species. 

6. The catalogued file, FUELS0L, exists and contains the fuel vapor mass 
fraction distribution determined previously by the VAPDIF code. 
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Then the following runstream is sufficient to execute the VAPDIF code when using the 
autoignition model. 

0ASG.AX EIGHT, D/0/TRK/300000 
0ASG.AX NINE, D/O/TRK/ 2 50000 
@USE 8, EIGHT 
@USE 9, NINE 

0ASG.T 10,D/0/TRK/2250000 
0ASG.AX ELEVEN, D/0/TRK/250000 
@USE 11, ELEVEN 
0ASG.T 12/D/0/TRK/750000 
0 ASG , T 13.D/0/TRK/2250000 
0ASG.AX SEVENTEEN., D/O/TRK/ 300000 
@USE 17 .SEVENTEEN 
0ASG.AX FUELS0L. , D/O/TRK/ 2 50000 
@USE 19.FUELS0L 

0ASG.AX TWENTYSEVEN. .D/O/TRK/300000 
0USE 27 , TWENTYSt'vEN 
0ASG.AX TWENTYEIGHT. .D/0/TRK/300000 
@USE 28 , TWE'’TYEIGHT 

@XQT . DIFFUSEMAP 
(INPUT CARDS) 

0FREE 8 
0FREE 9 


0FREE 10 
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0FREE 11 
0FREE 12 
0FREE 13 
0FREE 17 

@FREE 19 
0FREE 27 


(dFREE 28 
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6.2 Input Format for VAPDIF Cod® 

The Input to the VAPDIF code is described on the Input data coding forms 
which follow. These coding forms are arranged with one form per Input data card. 
In general the input data is read as follows: 

Card 1 Title Card 

Card 2 Option Caro 

Card 3 Print Option Card 

Card 4 Miscellaneous Data Card 

Card 5 Autoignition Model Constants (4) Cards 
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6.3 Output Description for VAPDIF Code 

The output on each page from the VAPDIF code is largely self-explanatory. 

A general description of the output by page is given below. 

Input Page 

The input page presents all of the input data including: options, computa- 

tional grid size controls, iteration parameters, and the axial coordinates of the 
corner point (J ■ 1, K « 1) . In addition the reference conditions for the ADD code 
viscous flow field are printed. 

Source Term Pages 

If IPRNT4*1, the source terms on the cvcssplane grid are printed These terms 
are arranged by (K row), (J column) where K is the index for the tangential (azimuthal) 
coordinate and J is the index for the normal (radial) coordinate. 

Solution Page 

If IPRNT1 = 0, the mass fractions on the crossplane grid are printed. 

Coordinate Grid Page 

For IPRNT2 or IPRNT3 t 0, the coordinates on the crossplane grid are printed. 
Iteration Page 

If IDBi - 1, the iteration history of the calculation is printed. This print- 
out includes the iteration number, residual, maximum concentration, and mass flow 
weighted average concentration for each residual. 

Flow Summary Page 

This page presents values of vapor flow rate and vapor fuel to air ratio as 
a function of axial location determined by integrating the vapor concentration over 
the crossplane grid at each axial station. The last column, labeled CTIL, is the 
overall vapor fuel to air ratio as determined by the PTRAK code and serves as a 
check of the accuracy of the VAPDIF calculation. 
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6. A Diagnostics for VAPDIF 

At the present time no diagnostics exist for the VAPDIF code. 
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6.5 Debug Options for VAPDIF Code 

If the options IDB1, IDB2, IDB3 ere set equal to unity, Intermediate results 
are printed for each iteration. 

Optlo ' Data Printed 

IDB1 Print maximum iteration residuals 

IDB2 Not used 

IDB3 Print coefficients of linearized equations 

The solution for each iteration at the point (J1DBG, K1DBG) on the calculation 
plane (crossplane grid) may also be printed. These terms are arranged by (K row), 
(J column) where K is the index for the tangential (azimuthal) coordinate and J is 
the index for the normal (radial) coordinate. Solutions for the points (J2DBG, 
K2DBG) , (J3DBG, K3DBG) , and (J4DBG, K4DBG) may also be displayed. 

It is recommended that IDB1 be set equal to unity. 
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6.6 Sample Input for VAPDIF Code 

The sample of Input to the VAPDIF code is based on the Swirl Tube Premixing 
Passage case described in Volume I, Section 7. Since almost all input data required 
to run the code is stored in data files generated by the ADD and PTRAK codes, 
little input is required by the VAPDIF code. On card 2 (line 2), it is noted that 
the number of iterations is limited to I0PT1-5. The initial flovfleld station is 
located at the first ADD code coordinate station (IADD“1) , the first calculation 
station after the initial plane is IBEGIN»2, and the calculation will terminate 
at station IEND = 5. The print option card (line 3) indicates that both the con- 
centrations and fuel source distributions will be printed at each station. The last 
card (line A) shows that the Schmidt number is unity. 
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7.0 GLOBAL STRUCTURE 


7.1 Interaction of ADD, PTRAK, VAPDIF Codes 

The Global Structure Flow Chart which describes the Interactions of the ADD, 

PTRAK, and VAPDIF codes Is shown on Fig. 7.1. The three coder are executed 
independently In the sequential order shown and output stored on data files. The 
data files (Units 8, 9, 11, 17, 19, 27, 28) are the only Interfaces between the 
codes. 

The ADD code calculates the computational coordinates (S,N) and physical 
cylindrical coordinates (r,z) and stores this data on Unit 9. This coordinate 
data is required by both the PTRAK and VAPDIF codes. The ADD code also calculates 
the viscous flow solution and stores the flow variables (P,T,p ,$,u E ) on Unit 8. 

These flow variables are required by both the PTRAK and VAPDIF. Certain control 
parameters, such as the number of streamlines KL and number of streamwise stations 
JL, which are required by the PTRAK and VAPDIF codes are also stored on Units 8 
and 9 and are not required inputs to these codes. This procedure reduces conflicts 
or ambiguities when executing the three codes. 

The PTRAK code calculates the droplet trajectories and degree of evaporation 
from the initial droplet conditions and known flow field. This solution is used 
to calculate the source terms Wf and Wf w which are stored on Unit 17 for use by the 
VAPDIF code. The PTRAK code also calculates preignition reactions and stores a 
source term Sj on Unit 27. As the droplets evaporate, they cool the surrounding 
air. This temperature drop AT is stored on Unit 28. 

The VAPDIF code can be executed in two modes. The first execution, with 
ISDG-0, calculates the three dimensional fuel concentration throughout the flow 
field using the fuel source distribution calculated from the PTRAK code. This 
solution is written on Unit 11. On the second execution, ISDG«1 or 2, VAPDIF 
calculates the concentration of critical species. The fuel concentration which was 
previously written on Unit 11 is read on Unit 19 and used to calculate the source 
term due co chemical reaction in the gas phase. The critical species concentration 
is then stored on Unit 11. 

It is noted that all control parameters are calculated and stored by the first 
code and are not required input to succeeding codes. This arrangement of codes produces 
considerable flexablllty in that any data which is calculated and permanently stored 
on data files need not be repeated. As an example: for any given premixing duct 
and airflow conditions, several different droplet ejectors may be examined using the 
PTRAK code without repeating the ADD code calculation. 
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7.2 Input/Output Data Files 

The output data files for the ADD; PTRAK, and VAPDIF codes are described 
on Tables 7.1, 7.2, 7.3 respectively. These tables show the Unit number, the 
names of the arrays stored on the file, the block (record) length in words, the 
number of blocks (records), and the subroutine generating the files. All unit 
numbers are set by parameters statements (in brackets). All arrays are single 
precision except AFF generated by the ADD code. 

In the operation of the ADD code, only files assigned to Units 8, 9, 11, 
and 22 are required. The other files are only used for special options in the 
ADD code which would not generally be used for a LPP duct calculation. Note also 
that only output files are shown under the computer code name. Input files for 
the PTRAK code are described on Table 7.1, and input files for the VAPDIF are 
described on Tables 7.1 and 7.2. The actual file assignments and 'unstreams 
for each code are described in Sections 4.1, 5.1, and 6.1. 

These data files are read/written by general 1/0 routines DR0UT in the ADD/81 
and \ APDIF codes or by BLKRED in the ADD/82 (Ref. 4) and PTRAK. These subroutines 
use :he UNIVAC Library I/O routine NTRAN. Both DR0UT and BLKRED are easily con- 
verted to ANSI standard F0RTRAN DEFINE FILE. 
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Table 7.1 

Table of File Assignments ADD Code 


UNIT 

ARRAY 

BLOCK LENGTH 

NO. 

SUBROUTINE 

NO. 

NAME 

(WDS) 

BLOCKS 

WRITING BLOCK 

8 (NDRUM) 

F(NEQ, 3, 1ST) 
FPARM(15) 

3015 

JL-2 

S0LVI 

9(JDRUM) 

JiTEP 

Q(19, 1ST) 

RHS(10) 

RMS(10) 

RTS (10) 

DSTEP 

QPARM(9) 

1941 

JL 

C00RST 

lO(CDRUM) 

FF(17 ,2 , 1ST) 

3400 

1 


ll(LDRUM) 

AFF(LNGTO)** 

6000 

5 

S0LVI 

12(LF0RC) 

F0RC 

780 

1 

F0RCE 

19(KP0IS) 

JSTEP 

Q( 19 , 1ST) 

RHS(IO) 

RMS (10) 

RTS(IO) 

DSTEP 

QPARM(9) 

1941 

JL 

C00RST 

22 (MDRUM) 

FIV(NEQ, 3, 1ST) 
FIPARM( 15) 

3015 

JL-2 

CAL IN V 

23(NFDRM) 

F ( 1ST4) 

104 

JL+NST-2 

P0ISCF 

24 (NPDRM) 

P( 1ST) 

100 

JL+NST-2 

P0IS 

2S(NGDRM) 

BLK(NGIST) 

400 

JL+NST-2 

P0ISCF 


**Double precision variables 
IS - 100 
NST “ 25 
1ST - 100 
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Table 7.2 

Table of Flic Aealgnmente PTRAK Coda 


UNIT 

ARRAY 

BLOCK LENGTH 

NO. 

SUBROUTINE 

NO. 

NAME 

(WDS) 

BLOCKS 

WRITING BLOCK 

12(KTDRUM) 

PART(200,13) 

2600 

JL 

SMTER 

17(MPRUM) 

AMASS( 1ST, 50) 

RDUM(IST,4) 

APARM(IO) 

5410 

JL 

DATAM 

18(NB0UN) 

RPART( 1ST, 2) 
PHPART(IST,2) 

400 

JL 

SMTER 

19( 1SDRU) 

SARRAY ( 8) 

8 

JL 

SMTER 

27(MDRUK1) 

AMASS (1ST, 50) 
RDUM( 1ST, 4) 
APARM(IO) 

5410 

JL 

DATAM 

28(MDRUM2) 

AMASS (1ST, 50) 
RDUM( 1ST ,4) 
APARM(IO) 

5410 

JL 

DATAM 


.11, 5 IS - 100 
KL i 1ST =■ 100 
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Table 7.3 

Table of File Assignments VAPD1F 


UNIT 

NO. 

ARRAY 

NAME 

BLOCK LENGTH 
(WDS) 

NO. 

BLOCKS 

SUBROUTINE 
WRITING BLOCK 

lO(IYDRM) 

Y (3 , 3, 1ST, 10) 

45,000 

JL 

APDC0R 

ll(IFDEM) 

FG( 1ST , 10) 

5,000 

JL 

S0LVE 

12 ( IXDRM) 

X( 3, 1ST, 10) 

15,000 

JL 

ADDC0R 

13( ITC0R) 

TC0R(3, 3, 1ST, 10) 

45,000 

JL 

ADDC0R 
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8.0 DETAIL DESCRIPTION OF PTRAK CODE 


8.1 Main Program/Flow Chart 

The flow chart for the main program PTRAK is shown on Fig. 8.1. It consists 
of three major tasks: 1) Read input data and set up initial conditions for the 

droplet (subroutines INPUT, INTIAL, 0UTPUT) ; 2) Calculate the droplet trajectory, 
evaporation rate, and production of critical species (subroutines TRACK, SUMRY) ; 
and 3) Calculate the source terms for use by the VAPDIF (subroutine DATAM) . 

In task 1, INPUT reads the input data according to the input format described 
in Section 4.2. INIT determines the initial droplet conditions (location, tem- 
perature, velocity, mass) which define the droplet classes by specifying number 
density distribution described in Section 4.7. 0UTPUT prints the input data and 
initial conditions. 

The integration of the droplet equations for one step is accomplished by 
subroutine TRACK using a fine coordinate grid interpolated from the course coordinate 
grid calculated by C00R. The outer D0 loop indexes the course grid where the J 
index corresponds to the Jth streamwise coordinate (block) stored on file 9. The 
inner D0 loop indexes -.he fine grid in KDS interpolated steps using the index JKDS. 

The integration takes place from JFIRST to JLAST in KDS*(JLAST-JFIRST) steps. 

At the present time there is no algorithm for varying the step size KDS. As the 
integration proceeds, the solution at each J ctl station for the droplet dependent 
variables is stored on file 12 and a summary of this data is stored on file 19. 

When the integration is complete, a summary of the solution is printed by SUMRY. 

The execution of the third major task depends on the input option I0PT7. 

When I0PT7>O, subroutine DATAM searches through the solution files (8,9,12) calculates 
the source terms for the production of fuel, the source terms for the production 
of critical species, and the drop in air temperature due to droplet evaporation. 

These results are stored for each mesh point on the course grid corresponding to the 
coordinate grid stored on file 9. 
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8.2 Description of Modules by Major Task 

A complete list and brief description of all the subroutines is given in 
Section 8.3. The flow chart for the main program PTRAK is shown on Fig. 8.1 
and an overall description is given in Section 8.1. The description of the 
module function? will follow this flow chart. Only major computational tasks 
will be describee Minor tasks and Input/Output are self explanatory and will 
not be described. 

Following through the flow chart Fig. 8.1, the first major task is performed 
by subroutine INTIAL which stores the crossection boundary at the initial station 
on unit NB0UND, reads the flow field variables at the initial station from unit 
NDRUM, and calculates the initial droplet conditions by calling subroutine FN0Z. 

FN0Z specifies the initial conditions (Eqs. 4.7.1), and specifies the number density 
in phase space defined by (Eq. 4.7.2). Generally the binomial distribution (Eq. 

4.7.3) calculated by FBIN is used. At the user's option, the Rosin-Rammier dis- 
tribution (Eq. 4.7.10) may be used to determine the distribution of droplet diameters. 

The second major task is integrating the droplet equations and calculating the 
droplet collision interactions. Subroutine TRACK performs this calculation for one 
time step. The flow chart for TRACK is given on Fig. 8.2 which shows the breakdown 
of the major task into eight smaller tasks. At the beginning of each station on the 
course grid (J station) IPR-1, the droplet data is read from unit KTDRUM and the 
flow field variables from unit NDRUM. For each step on the fine grid (JKDS>1) and 
IPR^O, and the flow field and coordinate variables are interpolated. 

The droplet class is defined by the index IJK calculated from the nested D0 
loop indices. Thus for each droplet, the droplet motion is calculated by PARPAT, 
a determination is made if the droplet hits a boundary by WALLRB, and depending 
on the options I0PT1 and IBREK, a determination is made on whether the droplet 
shatters by BREKUP, and whether the droplet is a member of the two largest classes 
by L0CGRD . When the solution for all the droplets at the new time step is determined, 
L0GC0L examines the two largest classes and applies the collision model. The 
solution is then printed by subroutine PRINT and stored by subroutine SMTER. 

The equations of droplet motion are contained in PARPAT which are integrated 
for one time step using the predictor corrector method described in Section 4.11. 

Table 8.1 shows a breakdown of the tasks performed by PARPAT. FINTP locates the 
droplet on the computational grid and interpolates the flow field variables at that 
point. FPR0P , GPR0P, and PPR0P calculates the thermodynamic and transport properties 
for the film, air, and, liquid respectively. DIFFUS calculates the mass diffusion 
coefficient (Eq. 4.4.16), P0LY is a polynomial interpolation, C0XCH is the Cox chart 
(see Section 4.5). DRAGF contains the drag coefficient correlation (Eq. 4.2.8) 
and NUSSET contains the Nusselt number correlations (Eqs. 4.4.1, 4.4.2) for heat 
and mass transfer. VAHR calculates the heat and mass transfer rates (Eqs. 4.3.4 and 
4.3.1). The equations of motion (Eqs. 4.2.1 through 4.2.6) are programmed directly 
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using DRAGF. The role of change of droplet temperature (Eq. 4.3.3) and droplet 
diameter (Eq. 4.3.9) is calculated by DTEMP and DDIAM respectively. When the 
predictor corrector iteration converges (Eq. 4.11.10) the time step is known 
and the calculation of critical species is made. X2INIT calculates the initial 
distribution of critical species in the film surrounding the droplet. X2INIT 
requires VA11R and PFILM to calculate the local partial pressure. XX2 integrates 
the rate equation (Eq. 8.3.1) using XRATE. X2MASS integrates over the film volume 
to obtain the mass of critical species (see Eq. 8.2.6). Finally the total heat 
input to the droplet (see Eq. 8.2.11) is calculated. 

The third task in TRACK (Fig. 8.2) is to examine the droplet to see if it 
hits a boundary using subroutine WALLRB. For purely elastic wall rebound, the 
velocity component normal to the wall is changed in sign. For periodic boundary 
conditions, a droplet leaving one boundary enters on the opposite boundary. A 
portion of the droplets may stay on the wall and form a liquid film. These con- 
ditions are calculated by subroutine B0UNCE. B0UNCE determines the fraction of 
droplets that remain on the wall using the models described in Section 4.10. Of 
the fraction that remains on the walls, some portion will evaporate depending on 
the wall temperature and partial pressure calculated from the Cox chart using 
subroutine C0XCHT as described in Section 4.10. 

The next step is to determine if a droplet shatters using subroutine BREKUP 
and the models discribed in Section 4.9. The droplet shattering model is applied 
to all classes. After shattering the droplets are counted with the neare- t class 
rather than forming new classes. However, the two largest classes in each com- 
putational grid must be recalculated using L0CGRD. 

Once the new droplets conditions are established for all classes (IJK loop 
completed), the droplet collision model L0GC0L can be applied. L0GC0L searches 
through the computational grid and determines if the two largest classes are in 
the same grid volume. If they -ire in the same grid volume, subroutine C0LLDE 
determines if a collision occurs, and the conditions after a collision using the 
model described in Section 4.8. A fraction of the droplets may coalesce and 
a fraction may rebound with elastic collisions. Of those that rebound, subroutine 
C0LLSN calculates the approach velocity along the line of centers and subroutine 
C0LDYN calculates the rebound velocity. Subroutine C0LLSN contains the velocity 
transformation matrix and its inverse to go from the computational coordinates to 
the collision coordinates. Following the collision, the new properties of the classes 
are calculated in C0LLDE. The solution is printed by PRINT and stored on files 
NBQUND, ISDRV,KTDRUM, by SMTER at each station of the course computational grid. 
This last step completes the calculation in TRACK. 

Returning to PTRAK, Fig. 8.1, the calculation continues until the outer D0 
loop is completed. An output summary is printed for all J stations on the course 
grid by subroutine SUMRY. Depending on the input option I0PT7, the source terms 
are calculated by DATAM. DATAM searches through all the solution files and calculates 
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• • 

the source terms Wf end for the vapor diffusion equation (Kq. 5.2.1) using 
Eqs. 4.7.14 and 4.7.15. iKie result Is stored on fils MDRUM. Zt then calculates 
the source term Sj for the critical species equation (Eq. 8.2.1) using Eq. (8.2.6) 
and stores the result on file MDRUM1 . Finally it calculates the temperature drop 
due to evaporlzation and stores the results on file MDRUM2. At this point the 
calculation In PTRAK Is complete. 
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Table 8.1 


Tasks in Subroutine PARPAT 


Subroutine 


Aux. Subroutines 


FINTP 


locate droplet on grid 

find local air flow conditions 


FPR0P 


NUSSET 


film properties (Eqs. 4.4.6 to 4.4.20) 
air properties 

liquid properties (Eqs. 4.4.21 to 4.4.28) 

drag coeff. (Eq. 4.2.8) 

Nusselt No. (Eq. 4.4.1, 4.4.2) 

determine heat transfer (Eq. 4.3.4) 
determine mass transfer (Eq. 4.3.1) 

Equations of motion (Eqs. 4.2.1 to 4.2.6) 

Rate of temperature change (Eq. 4.3.3) 
Rate of diameter change (Eq. 4.3.9) 


DIFFUS 

P0LY 

POLY, C0XCH 


X2INIT 


X2MASS 


Initialize critical species 
Integrate rate equation (Eq. 8.3.1) 
Integrate mass (bracket Eq. 8.2.6) 


VAHR, PFILM 
PFILM, XRATE 


Integrate bracket (Eq. 8.2.11) 
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ALJP0T 

BLKDAT 

BLKRED 

B0UNCE 

BREKUP 

C0LDYN 

COLLDE 

C0LLSN 

C00R 

C0XCHT 

DATAM 

DDIAM 

DIFFUS 

DISTEM 

DRAGF 

DTEMP 

FB1N 

FINTP 

FN0Z 

FPR0P 

GAMFUN 

GPR0P 


8.3 List of Subroutines in PTRAK 

Evaluate Leonard Jones Potential collision integral 
Block data 

Read/Write data on files 
Wall rebound model 
Droplet Shattering model 

Calculate droplet velocity after collision 
Droplet Collision model 

Calculate velocity transformation matrix 
Read coordinate data 

Determine vapor pressure from Cor chart 
Store fuel sources for VAPDIF 
Calculate rate of change in droplet diameter 
Calculate mass diffusivity 
Determine distillation temperature 
Calculate droplet drag 

Calculate rate of change of droplet temperature 

Calculate binomial distribution 

Interpolate flow field of air 

Calculate initial conditions for droplets 

Determine fluid properties of film 

Calculate Gamma function 

Determine gas properties for droplet 
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INPUT 

Read Input data 

INTIAL 

Initialisation of partlcla clasaaa 

L0CGRD 

Determine two largest dasaes 

L0GC0L 

Logic for Droplet collision model 

L00K 

Interpolation for table lookup 

NUSSET 

Calculate Nusselt number 

0UTPUT 

Write Initial input conditions 

PARPAT 

Integrate particle path equations 

PFILM 

Calculate local film properties 

P0LY 

Perform polynomial curve fit 

PPR0P 

Determine fluid properties of droplets 

PRINT 

Print class properties at each station 

PTRAK 

Main program 

SMTER 

Calculate summary terms 

SUMRY 

Write summary output 

TRACK 

Track droplet particle classes 

VAHR 

Calculate heat and mass transfer rates 

WALLRB 

Logic for wall rebounds 

XRATE 

Calculate net production of X 2 

XX2 

Calculate X 2 for one time step 

X2INIT 

Initialize X 2 mole concentration 

X2MASS 

Calculate mass of X 2 in droplet 
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8. A List of C0MM0N BL0CK Variables 


The C0MM0N BL0CK variables are grouped into labeled C0MM0N BL0CKS. 
Alphabetical listing if by labeled C0MM0N BL0CK name. The C0MM0N BL0CK name 

is given and a genet al description of the variables in the group. Following this 

is a detailed list of variables in the order in which they appear in the C0MM0N 

BL0CK. C0MM0N BL0CK variables for the ADD code are given in Ref. A. 

C0MM0N/AC0NS/ 

ADD 

Code Variables (Ref. A) 

COMMON /APART/ 


Droplet Input Conditions 

RNZCLI(N,1) 

r 

Injector radius (cm) 

RNZCLI(N,2) 

z 

Injector axial location (cm) 

RNZCLI(N, 3) 

♦ 

Injector circumferential location (rad.) 

RNZI(N) 


Radial displacement of r (cm) 

ALPHI(N) 

a 

Angle between normal and streamwise (deg) 

BETAl(N) 

6 

Angle between tangential and streamwise (deg) 

VMF.ANI(N) 

V 

Magnitude of velocity (K/sec) 

DPRTI(N) 

°L 

Mean diameter of droplet (micron) 

DDPRTI(N) 

6D l 

Variance of droplet diameter (micron) 

NN0Z 


Number of ejectors 

DALP(N) 

5u 

Variance in a (deg) 

DBET(N) 

66 

Variance in 6 (deg) 

IL0C 


Number of ejectors 

IVS 


Number of velocity classes 

IPHI 


Number of a classes 

ITHE 


Number of 6 classes 

IDIA 


Number of droplet diameter classes 

DVMI(N) 

6V 

Variance of Mean Velocity 


auBiililtti 
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COMMON / BONVAR/ Crowe tion Boundary Variables 


C1,C2 


Wall rebound probability conatanta 

RTIP 

r H 

Tip radius (cat) 

RHVB 

r T 

Hub radius (cm) 

LPHI 


Number of $ grid point 

DPHI 

A* 

Increment in 4 (rad) 

LPART 


Index in $ 

PHIL 

♦l 

Left $ boundary (rad) 

PHIR 

♦r 

Right 4 boundary (rad) 

RPART(L.l) 

r t (*> 

Tip boundary (cm) 

RPART (1,2) 

r H<») 

Hup boundary (cm) 

PHPART (K,l) 

♦L^ r ) 

Left & boundary (rad) 

PHPART (K,2) 

♦ r (r) 

Right & boundary (rad) 

COMMON/BPART/ 

Droplet 

Input Thermodynamic Properties 

RH0LPI(1) 

PL 

Droplet density (gm/cm^) 

TEMPI (I) 

t L 

Droplet temperature (*K) 

DLP(I) 


Droplet diameter (micron) 

KPART(N) 


Radial index for N** 1 droplet 

PERL0C(I) 

P L 

Percentage of fuel 

JPLL 
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C OMMON /CARRAY/ Thermodynamic Polynomial Coeff. 


C0EFP(1,I) 

A I 

Droplet density (gm/cm^) 

C0EFr(2,I) 

A I 

Droplet heat capacity (cal/cm/'K/sec) 

C0EFP(3,I) 

A I 

Droplet viscosity (gm/cm/sec) 

C0EFP(4 , I) 

A I 

Gas heat capacity (cal/gm/*K) 

C0EFP(5 , I) 

A I 

Gas conductivity (cal/cm/ # K/sec) 

C0EFP(6,1) 

A I 

Gas Viscosity (gm/cm/sec) 

C0EFP(7,I) 

A I 

Air heat capacity (cal/gm/*iC) 

COEFP(8, I) 

A I 

Air conductivity (cal/ cm/ *K/ sec) 

C0EFP (9,1) 

A I 

Air viscosity (gm/cm/sec) 

F(T) 

4 

-L A i 

i 1 - 1 


1-1 


COMMON CCRATE Parameters for Critical Species 


AB, AF 

A b* A f \ 

EB, EF 

E b* E f 1 

ALPB, ALrF 

a b' a f \ Constants in rate equation 

BETB, BETF 

8b# 8f 1 (Eq. 8.3.1) 

CAMB, GAMF 

y B* >B 1 

CARBN 

6 ' 

M0LE1 ,M0LE2 

,M 2 Molecular weights 

BRAD 

( r < j + b)/ r d Film thickness ratio 

Y02 

Yq 2 Mole fraction of O 2 

DELM2 

(X 2 /Xi) 0 Initial ratio 

LRLP 

Number of points in film 
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DELTB , DELTF 


Constants for rsts equation 

ANUB.ANUF 

% v f 

(Eq. 8.3.11) 


Cjjfggt /CCX21/ Mole Concentration Critical Species 


X21 

x 2 (t) 

Mole concentration (nole/cn^) 

X22(I) 

X 2 (t+At) 

3 

Mole concentration (mole/cm ) 

DX2DO(I) 

(dX 2 /dt) 

Production of X 7 (aole/cm^/sec) 
t+A L 

DX2DT(I) 

(dX 2 /dt) 

t Production of X 2 (mole/ cm-*/ sec) 

C0MM0N /COLLVA/ 

Index for 

Largest Class 

KGE0M(K,L,1) 


Index for largest class 

KGE0M(K,L,2) 


Index for second largest class 

Where K = r grid index 


L = $ grid indei 


C0MM0N /C0NST/ 

ADD Code 

(Ref. 4) 

COMMON /CONVER/ 

Constants 


GC 

g 

Force equivalent of mass (1.01325x10^ dyn/cm/atm) 

JC0NST 

J 

Energy equivalent of work (41.311cm mole/ atm/ °K) 

R0C0N 

R 

Universal Gas Constant (82.087 cm-*atm/mole/°K) 

R0C0NP 

R 

Universal Gas Constant (1.98717 cal/mol/°K) 

C0MM0N /CORE/ 

ADD Code 

(Ref. 4) 

C0MM0N /CORE 2/ 

ADD Code 

(Ref. 4) 

C0MM0N /COXCUR/ 

Cox chart data 

IC0X1 


Number of input points 

PICH 

P 1 

Vapor pressure (atm) 
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P2CH 

P 2 

Vapor pressure (atm) 

C0XCAR(1 , I) 

n 

Carbon number 

C0XCAR(2 , 1) 

T l 

(T(;,P L ) ( °K) 

C0XCAR( 3, I) 

T 2 

T( 4 * P 2 ) (°K) 

C0MM0N /CPART/ 

Dependent droplet variables 

PART( 1,1) 

U s 

Streamwise velocity (m/sec) 

PART( 1,2) 

U * 

Tangential velocity (m/sec) 

PART( I, 3) 

U n 

Normal velocity (m/sec) 

PART( I, 4) 

n 

Normal coordinate 

PARK 1, 5) 


Circumferential coordinate 

PART( I, b) 

d l 

Droplet diameter (micron) 

PART( I, 7) 

t l 

Droplet temperature (°K) 

PARK 1, 8) 

r 

Droplet radius (cm) 

PART (I, 9) 

z 

Droplet axial distance (cm) 

PART( I, 10) 

t 

Droplet time (sec) 

PART( 1,11) 

f 

Probability 

PARK 1,12) 

m 2 

Mass of critical species (gm) 

PART (I, 13) 

Q 

Total heat load (cal) 

C0MM0N /C_URV_K/__ 

D is t ilation 

curve 

IDSTIL 


No. of data points 

DISCUR(l.l) 

P e 

Percent evaporated 

DISCUR( 2 , I ) 

T b 

Temperature (°K) 
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C0MM0N /DREDl/ ADD code 

(Ref. 

4) 

C0JW0N /DSKL0C/ 

ADD code 

(Ref 

. 4) 

C0MM0N /DUC0UT/ 

ADD code 

(Ref , 



C0MM0N 

/DUCTIN/ 

ADD code (Ref. 4) 

C0MM0N 

/FILVAR/ 

Film thermodynamic properties 

MF 




Film molecular weight 

RH0F 


p m 


Film density (gm/cm^) 

CPF 


C 

pm 


Film heat capacity (cal/gm/°K) 

KF 


k 

m 


Film thermal conductivity (cal/cm/°K/sec) 

MUF 


l ‘m 


Film viscosity (gm/cm/sec) 

DMASSF 


m 


2 

Film mass diffusivity (cm /sec) 

CFF 


C f 


Fuel mass fraction 

CGF 


C 

a 


Air mass fraction 

SIGP6F 


0 

m 


O 

Film force constant (A) 

OMEGF 


ft* 


Collision integral 

YINF 




Air mole fraction 

YGF 


X - Y m 


Average film air mole fraction 

YPF 


Y m 


Average film fuel mole fraction 

YPSF 


Y f.s 


Fuel mole fraction at droplet surface 

CKPGF 


(e/<) 

m 


Force constant film (°K) 

TSTARF 


* 

T 


Dimensionless temperature 

YPINF 


V 


Mole fraction of air at film edge 

MINF 


M 

a 


Molecular weight air 
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C0MM0N /FLAGS/ ADD code (Ref. 4) 


C0MM0N /FUNC/ 

ADD code (Ref. 

JO 


C0MM0N / GASVAR/ 

Gas thermodynamic properties 


MG 

M 

a 

Air molecular weight 


PCG 

P 

ca 

Air critical pressure (atm) 


TCG 

T 

ca 

Air critical temperature ( 8 K) 


SIGG 

G 

a 

O 

Air force constant (A) 


CKG 

(e/<) 

a 

Air force constant (°K) 


CPG 

C 

pa 

Air heat capacity (cal/gm/°K) 


KG 

k 

a 

Air thermal conductivity (cal/cm/ °K/sec) 


MUG 

u 

a 

Air viscosity (gm/cm/sec) 


C0MM0N /INTINP/ 


Input Flags 


KL 


No. of streamlines 


JL 


No. of streamwise stations 


KDS 


No. steps/course grid 


KLL 


Not used 


NB 


Not used 


ISHAPE 


Not used 


I0PTN 


Input option see Section 5.2 


N = 1,17 




C0MM0N /MASSD/ 

Source terms 

for VAPDIF code 


AMASS(K,L) 

0 

w ,S ,AT, 

Source terms for VAPDIF 


RDUM(K.ISIDE) 

0 

w. , 

fw 

Wall source terms for VAPDIF 



M2-9 15362-40 


APARM(l) 


» 

Circumferential step else (red) 

APARM(2) 

LPHT 

• 

No. circumferential stepa 

APARM(3) 

U 

fv 

• 

Weight flow fuel (go/sec) 

APARM(4) 

w_ 

fa 

» 

Weight flow air . (gu/ sec) 

APARM(5) 

JLAST 

% 

Number of blocks 

APARM(6) 

I$PT3 

9 

mj'ut option 

APARM(7) 

JC$UN 

> 

Block routine 

COMMON /MISCEL/ 

Miscellaneous variat Les 

VEL 

V 


Air velocity (m/sec) 

TEMG 

T 

a 


Air temperature (°K) 

PRESG 

P 

a 


Air pressure (atm) 

RH0G 

P 

a 


3 

Air density (gm/cm ) 

C0MM0N /NUMDEN/ 

Droplet Summary 

FMASS 

0 

h l 


Fuel flow rate (gm/sec) 

FSUM 

n t 


No. fuel droplets/sec (1/sec) 

SMDD 



SMD denominator 

SMDN 



SMD numerator 

SMD 

SMD 


Sauter mean diameter (microns) 

QF 

0 

W L 


Liquid fuel flow rate (gm/sec) 

QFO 

O 

w 

LO 


Initial liquid fuel flow rate (gm/sec) 

SARRAY(l) 

z 


Axial location (cm) 

SARRAY ( 2) 

SMD 


Sauter mean diameter (micron) 

SARRAY ( 3) 



% fuel liquid 
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Ilk 


SARRAY(4) 


% fuel at wall that Is liquid 

fl 

SARRAY(5) 


% fuel at wall that Is vapor 

1] 

SARRAY(6) 


7 . fuel in vapor state 

{ I 

SARRAY(7) 


X fuel evaporated 

0 

SARRAY ( 8) 


Flow rate of vaporized fuel 

n 

SARRAY(9) 


Fuel/Air ratio 

M 

n I 

C0MM0N /NUSVAR/ 

Nusselt 

numbei variables 

ill 

\ 

NUH 

Nu 

h 

Nusselt number heat transfer 

ill 

■ 

NUM 

Nu 

m 

Nusselt number mass transfer 

n! 

REN 

Re „ 

Reynolds number of droplet 

* > i 

i 

PRN 

Pr 

Prandtle number 

i! 

ft* 

SCN 

Sc 

Schmidt number 

j ! 

COMMON /PARVAR/ 

Droplet 

thermodynamic properties 

, : j 

MP 

\ 

Molecular weight fuel 

! 1 
. i 

TCP 

T 

c 

Critical temperature of fuel (°K) 

l 

t ‘ 

4 1 

PCP 

P 

c 

Critical pressure of fuel (atm) 

’ T 

SIGP 

°L 

O 

Force constant (A) 

“ ! 
_ | 

CKP 

u/ *\ 

Force constant (°K) 

0 

LAMBP 

\ 

Heat oi vaporization (cal/gm) 

1 

TT : 

H i 
_ _ 

SP 

S 

Surface tension (dyne/cm) 

] 

„ i 

PVP 

P f,s 

Vapor pressure fuel (atm) 

T 1 

IMP 

T 

m 

Film temperature (°K) 

r? 

RH0LP 

P L 

3 

Fuel liquid density (gm/cm ) 

ii 




i ' 

CPLP 

C PL 

Fuel liquid heat capacity (cal/gm/°K) 

« y 

; 
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MULP 
CPVP 
KVP 
MU VP 


p Fuel liquid viscosity (gm/cm/sec) 

L 

C Fuel vapor heat capacity (cal/gm/°K) 

k^ Fuel vapor thermal conductivity (cal/cm/°K/sec) 

U v Fuel vapor viscosity (gm/cm/sec) 


C0MM0N /QM0N/ 


Q0LD(I) 


0 

q si 

Heat transfer rate (cal/sec) 


DMDTO(I) 


(dm2/dt) i 

Mass rate of change of critical species 


C0MM0N 

/REALIN/ 

ADD code 

(Ref. 4) 


C0MM0N 

/SVARB/ 

ADD code 

(Ref. 4) 
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A list of 
alphabetically 

ALJP0T 

TS 

ALJP0T 

BLKDAT 

BLKRED 

B0UNCE 

IJK 

IB0UN 

SP 

TEMMIN 

TEMMAX 

ND 

PRB 

WFC0LL 

WFEVDP 

WFEVAT 

WFLIQT 

BREKUP 

IJK 

IPREK 

A, B 

AK1.AK2 


8.5 List of Local Variables 

selected local variables appears below. These variables are 
listed by subroutine name. 

T* , Dimensionless temperature 

* 

ft , Collision integral 

ADD code (Ref. A) 

ADD code (Ref. 4) 

Class index 
Wall index 

S Surface tension (dyne/ cm) 

Minimum droplet temperature (°K) 

Maximum droplet temperature (°K) 

Number of droplets (1/sec) 

Probability of wall rebound 

2 

Amount of fuel/area (gm/cm ) 

2 

Amount of fuel evaporated/area (gra/cm ) 
Total fuel evaporated (gm) 

Total fuel on wall (gm) 

Droplet class index 

Flag indicates shattering 

Constant Wolfe, Anderson shattering model 

Aerodynamic drag constants 




D 

0 

ii 

0 

11 

n 

ii 





-4 


? % 


»• 



-- * 


;i 

1 5 

• 

*7 





•** 


ll 
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CD 

C D 

Drag coefficient C n *l 
D 

DK3 


Shattering Model constants DK3-136 

FK1.FK2 


Friction drag constants FKl-2, FK2-4 

RH0LP 

P L 

3 

Density of droplet (gm/cm ) 

SP 

s 

Surface tension (dyn/cm) 

VEL 

|V| 

Absolute velocity (m/sec) 

DSHAT 


Droplet diameter after shattering (micron) 

ND 

n t 

Total number density after shattering (l/sec) 

TEST1 


Time of aerodynamic breakup (sec) 

TEST2 


Time of friction breakup (sec) 

C0LDYN 



IP1.IP2 


Index for colliding classes 

VEL11 , VEL12 , VEL13 


Velocity after collision (Vs,Vn,V$) (m/sec) 

VEL21.VEL22.V23 


Velocity after collision (Vs,Vn,V$) (m/sec) 

AMASS1 ,AMASS2 


Droplet mass before collision (gm) 

E,EI 

E,E -1 

Velocity transformation matrix 

VC0L1 , VC0L2 

-> -*■ 

V V 

V K 

Velocity in collision coordinates (m/sec) 

UC1.UC2 

V V 

CJ1* CK1 

Velocity before collision (m/sec) 

VCL.VC2 

V V 

CJ2* CK2 

Velocity after collision (m/Bec) 

C0LLDE 



K, J 


Droplet index 

RC 

R c 

Distance between droplets (cm) 

XC,YC,ZC 

* 

Relative coordinates (cm) 
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COLLSN 



IP1 , 1P2 


Droplet index 

XC,YC,Zc 

X ,Y ,Z 
c c c 

Relative cartesian coordinates (cm) 

A(l, J) 


Transform matrix streamline to cylindrical 

B( 1 , J) 


Transform matrix cylindrical to cartesian 

C(1,J) 


Transform matrix streamline to cartesian 

D(I,.I) 


Transform matrix cartesian to droplet 

E(I,J) 

E 

Transform matrix streamline to droplet 

EI(1,J) 

E‘ l 

Transform matrix droplet to streamline 

C00R 

ADD coda 

(Ref. 4) 

C0XCHT 



PRESS 

P 

V 

Vapor pressure (atm) 

BC0EF 

8 

n 

Cox chart constants 

TEM1 

T 

L 

Droplet temperature 

PERCT 

P 

c 

Percent evaporated 

DATAM 



I0P 


Module switch 

JC0UN 


Block counter 

DDIAM 



WD0TSF 

O 

W L 

Vaporization rate (gm/sec) 

DRDT 

d(p L )/dt 

3 

Rate of density change (gm/cm /sec) 

DLP 

d l 

Droplet diameter (micron) 

DDIAM 

d(D )/dt 

L 

Rate of diameter change (micron/sec) 
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DIFFUS 




TEMP 

T 

Temperature (°K) 


PINF 

P 

a 

Air pressure (atm) 


DIFFUS 


Mass diffusion coefficient 

(cm^/sec) 

DISTEM 




PER 

P 

e 

Percent liquid evaporated 


DISTEM 

t d 

Dlstallatlon temperature 


DRAGF 




RH0 

P L 

3 

Droplet density (gm/cm ) 


DLP 

°L 

Droplet diameter (micron) 


TINF 

T a 

Air temperature (°K) 


PINF 

P 

a 

Air pressure (atm) 


VINF 

v i 

Air velocity (m/sec) 


VP ART 

v i 

Droplet velocity (m/sec) 


VEL 

AV 

Relative velocity (m/sec) 


CD 

C D 

Drag coefficient 


DTEMP 




QNET 

o e 

w 

Net heat flux (cal/sec) 


DLP 

d l 

Droplet diameter (micron) 


DTEMP 

dT/dt 

Rate of temperature change 

(°K/ sec) 

FBIN 




IXX 


IXX th class 


JXX 


Total number of classes 


FBIN 

f(I.IL) 

Probability 
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FINTP 


ANE 

n 

Normal coordinate 

W 


Metric coefficient 

VS, VP 

V ,V 

8 4 ) 

Air velocity components (m/sec) 

R0 

p 

a 

3 

Air density (gm/cm ) 

TP 

T 

a 

Air temperature (°K) 

PS 

P 

a 

Air pressure (atm) 

RB,ZB 

r,z 

Coordinates of droplet (cm) 

RSB.RNB 


Direction cosines of coordinates 

VSB.VNB 


Curvatures of coordinates 

KP 


Streamline index for droplet 

M 


Index for solution variables 

FN0Z 



ANP 

n 

Normal coordinate 

HE IT 

H 

Duct height (cm) 

FL , FV , FT , FP , FD 


Probability functions 


f 

0 


pF 


Circumferential location (rad) 

RP ,ZP 

r,z 

Droplet coordinates (cm) 

VSP,VNP,VPP 

V ,v ,v 

s n p 

Droplet velocity (m/sec) 

FPR0P 

PINP 

P 

a 

Air pressure (atm) 

GAMFUN 

ARC 

X 

Argument 

GAMFUN 

r(x) 

Gamma Function 
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GPR0P 

CPC 


c 

p* 

Air heat capacity (cal/gm/*K) 

KG 


k 

a 

Air thermal conductivity (cal/ca/*K/aac) 

MUG 


w a 

Air viscosity (ga/cm/sac) 

INPUT 

See 

Section 5.2 


INTIAL 

See 

C0MM0N BLOCK variables 

LOCGRD 

I.JK 



Index for iroplet class 

L0GC0L 




X1,Y1,Z1 


W z j 

Cartesian coordinates droplet (cm) 

X2,Y2,Z2 


VY z c 

Cartesian coordinates K^ droplet (cm) 

XC,YC,ZC 


X e*V Z « 

Relative cartesian coordinates (cm) 

RC 


R 

c 

Distance along line of center (cm) 

L00K 




X(I) 


x i 

Table of abscissa values 

Y(I) 


Y i 

Table of ordinate values 

XIN 

• 

X 

Input abscissa 

Y0UT 


Y 

Output ordinate 

KK 



Default flag 

NUSSET 

See C0MM0N /NUSVAR/ 


0UTPUT 

Self explanatory 


PARPAT 




VS1.VN1,' 

VP1 

V * v ’ V * 
s n ♦ 

Droplet velocity components at t (m/sec) 

AN1 


n 

Normal coordinate at t 

PI 


* 

Circumferential location at t (rad) 

DLP1 


d l 

Droplet diameter at t (micron) 
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TLP1 

t l 

Droplet temperature at t (*K) 

VS2,VN2,lfP2 

V V V 
s n 4> 

Droplet velocity components t+dt (m/sec) 

AN 2 

n 

Normal coordinate at t+dt 

P2 

* 

Circumferential location t+dt (rad) 

DLP2 

d l 

Droplet diameter at t+dt (micron) 

TLP2 

t l 

Droplet temperature at t+dt (®K) 

RB2.ZB2 

r,z 

Droplet coordinates at t+dt (cm) 

RH0LPO 

P LO 

3 

Initial droplet density (gm/cm ) 

DLPO 

°LO 

Initial droplet diameter (micron) 

KP 

KP 

Streamline index 

OT2 

AC 

Time Increment 

QS1 

• 

q 

Droplet heat transfer rate at t (cal/sec) 

X2M2 

m 2 

Mass of critical species at t+dt (gm) 

QHEAT1 

Q 1 

Total heat added to droplet at t (cal 

QHEAT2 

Q 2 

Total heat added to droplet at t+dt (cal) 

DM1 

d(m2)/dt 

Rate of change of mass at t (gm/sec) 

PFILM 



RL 


Droplet film radius (micron) 

RLP 

f L 

Droplet radius (micron) 

TLP 

t l 

Droplet temperature (*K) 

PINK 

P 

a 

Air pressure (atm) 

PF 

p f 

Vapor pressure fuel (atm) 

IF 

‘f 

Temperature of fuel (*K) 

XF 

x r 

Mole concentration of fuel (mole/cm^) 
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P0LY 


I 


Property Index 

TEM 

r 

Temperature *K 

PPR0P 



T1NK 

T 

a 

Air temperature (*K) 

TLP 

\ 

Droplet temperature ( # K) 

Dll’ 

\ 

Droplet diameter (micron) 

KH0LPO 

°L0 

3 

Initial droplet density (Rm/cm ) 

DIMA) 

D 

LO 

Initial droplet diameter (micron) 

PRINT 

Self 

explanatory 

PTRAK 



KMC 


Conversion factor (cm/ft) 

KHUB.ZHDB 

V*H 

ID wall coordinates (cm) 

RT1P.ZTIP 

V*T 

OD wall coordinates (cm) 

SMTKR 



WA1 R 

0 

w 

a 

Weight flow air (Rra/uec) 

WFKVAT 


Weight flow gus (gra/sec) 

WFLIQT 

0 

\ 

Weight flow liquid droplets (gm/sec 

QK 

o 

u 

LO 

Initial fuel flow (gm/sec) 

z 

z 

Axial location (cm) 


SUMKY Self explanatory 
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TRACK 



1D1A 


No. diameter classes 

IL0C 


No. location classes 

I PHI 


No. normal angle classes 

1THE 


No. asmutial angle classes 

IVS 


No. velocity classes 

1JK 


Class index 

VAHR 



QNET 

O 

q -w X 
s s 

Net heat transfer rate (cal/sec) 

WD0TSF 

© 

W L 

Vaporization rate (gm/sec) 

TLP 

t l 

Droplet temperature (°K) 

DLP 

d l 

Droplet diameter (micron) 

T1NF 

T 

a 

Air temperature (°K) 

P1NF 

P 

a 

Air pressure (atm) 

WALLRB 



1JK 


Droplet class index 

XRATE 




T 

XI 

X2 

X02 


T 

X. 


02 


Temperature (°K) 

3 

Mole concentration species 1 (mole/cm ) 

3 

Mole concentration species 2 (mole/cm ) 
Mole concentration oxygen (mole/cm^) 
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XX2 



DLP1 

d l 

Droplet diameter at t (microns) 

DT2 

At 

Time step (sec) 

TLP 

t l 

Droplet temperature (°K) 

PINF 

P 

a 

Air pressure (atm) 

TINF 

T 

a 

Air temperature (°K) 

X2INIT 



DLP 

d l 

Droplet diameter (micron) 

TLP 

t l 

Droplet temperature (°K) 

PINF 

p 

a 

Air pressure (atm) 

TINF 

T 

a 

Air temperature (®K) 

X2M2 

M 2 

Mass critical species (gm) 

QIT0T 

9 l 

Total heat added to droplet (cal) 

X2MASS 



BRAD 

<r d +b)/r d 

Film thickness ratio 

DLP1,DLP2 

d l 

Droplet diameter at t Rt+dt (micron) 

Dn2DT 

dM^/dt 

Rate of mass change critical species (gm/sec) 

DT2 

dt 

Time step sec 
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9.0 DETAIL DESCRIPTION OF VAPDIF CODE 
9.1 Main Program/Flow Chart 

The flow chart for the main program for the VAPDIF code is shown on Fig. 9.1. 

It consists of four major tasks: 1)~ Read and print the input data; 2) Calculate 

the three dimensional coordinates; 3) Calculate or restart initial conditions; and 
4) Solve the diffusion equation. The run stream for execution of the code is given 
in Section 6.1. This run stream Includes both input and output data files. A 
table of output file assignments is given on Table 7.3 and a table of unit numbers 

is given by Table 9.1. Table 9.1 shows both input and output files. N0P is the 

unit index used by the I/O routines. For input files the code generating the data 
is given by the code name in brackets. For output files the subroutine generating 
the data is given without brackets. 
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9.2 Description of Modules by Major Task 

A complete list and brief description of all the subroutines is given in 
Section 9.3. The flow chart for the main program for VAPDIF is shown in 
Fig. 9.1 with an overall description given in Section 9.1. The description of 
the module functions will follow this flow chart. Only major tasks will be described. 
Input/Output routines are self explanatory and will not be described. Of the task 
shown on Fig. 9.1 only tasks two and four need be described. 

Subroutine ADDC0R generates the three dimensional coordinate system used by 
the VAPDIF. Coordinates generated by the ADD code and stored on file 9, are 
read block by block. Each block corresponds to on streamwise station. ADDC0R 
then calculates the coordinates, metrics, and arc length distance and stores the 
result on unit 10. It then calculates the three dimensional cartesian coordinates 
for each mesh point and stores the results on unit 12. Finally it calculates the 
transformation matrix from the computational coordinates to the cartesian coordinates 
and stores the result on unit 13. Data on units 12 and 13 are not necessary for 
the solution but are useful for plotting results in physical space. 

Subroutine S0LVE solves the diffusion equation (Eq. 5.2.1) or the critical 
species equation (Eq. 8.2.1) depending on the input option ISDG. The main program 
sets the first and last station. A flow chart for subroutine S0LVE is shown on 
Fig. 9.2. The first task is to align the data blocks on units 8, 17 and 28 with 
the absolute coordinate location of the coordinate grid stored on unit 9. The I D 0 
loop steps off streamwise stations from IL0W to INUM. Next flow field data calculated 
by the ADD code and stored on unit 8 is read and source terms stored on units 17 and 
28 are read. The V D0 loop is an iteration loop for the point relaxation algorithm 
described in Section 5.1, and the K and J D0 loops sweep the entire crossplane grid 
including the boundary points. The coordinates stored on unit 10 are read by RDINB. 
The coefficients of the differential equations are calculated by C0EFFI. For ISDG>0j 
C0EFSG calculates the source term for the critical species equation (Eq. 8.2.1). 

On the boundaries, PBNDC applies normal derivative boundary conditions, and PERBC 
applies periodic boundary conditions depending on the input options. P0ISSN applies 
the difference operators (see Section 5.1) and solves for the guess. With the 
completion of J, K grid sweep, a convergence check is made. If the solution con- 
verges, the solution is printed by PRTDSK and stores on unit 11 by WR0UT. The 
calculation then moves to the next station until the I D0 loop is completed. The 
calculation then returns to the main program. 
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Table 9.1 

Table of Unit Numbers 


N0P 

NAME 

UNIT 

SUBROUTINE 

(CODE) 

DATA 

1 

IYDRM 

10 

ADDC0R 

3D coordinates 

2 

IFDRM 

11 

S0LVE 

Output file 

3 

IXDRM 

12 

ADDC0R 

3D Castes lan coordinates 

4 

ITDRM 

13 

ADDC0R 

3D transformation matrix 

5 

INDRM 

17 

(VAPDIF) 

Mass fraction of fuel 

6 

IADDRM 

8 

(ADD) 

Air flow field 

7 

IC0DRM 

9 

(ADD) 

Axisymmetric coordinates 

8 

IFGDRM 

19 

(PTRAK) 

Summary droplet data 

9 

28 

28 

(PTlwJC) 

Critical species source 
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9.3 List of Subroutines in VAPDIF 


ADDC0R 

Calculate 3D coordinates 

AINTG 

Area integral of prescribed function 

C0EFFI 

Calculate coefficients of P.D.E. 

C0EFSG 

Calculate source for critical species 

DERBDY 

Calculate derivatives on boundaries 

DERINT 

Calculate interior derivatives 

DERIV 

Calculate derivatives 

DFC0R 

Correct fuel concentrate 

DIFFUSE 

Main program 

DM2DT 

Critical species rate equation 

DR0UT 

Input/Output module 

FL0WIN 

Initial conditions module 

PBNDC 

Neuman boundary conditions 

PERBC 

Periodic boundary conditions 

PERIOD 

Read boundary values 

P01SSN 

Solve Poisson's equation 

PRINPT 

Print module 

PRUN 

Metric conversion 

PRTDSK 

Print module 

RDINA 

Read contiguous block data 

RDINB 

Read overlapping block data 

RESTRT 

Restart module 
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RH0V 

S0LVE 

S0URCE 

SUMRY 

WR0UT 

WR0UTB 

WTFL0W 


Integrand for flow Integration 

PDE solution algorithm 

Integrand for area integral 

Print summary data 

Write contiguous block data 

Write overlapping block data 

Integrand of Mass flow area integral 
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9.4 List of C0MM0N BL0CK Variables 

The C0MM0N BLOCK variables are grouped into labeled C0MM0N BL0CKS. Alphabetical 
listing is by labeled C0MM0N BL0CK name. The C0MM0N BL0CK name is given and a 
general description of the variables in the group. Following this is a detailed 
list of variables in the order in which they appear in the C0MM0N BL0CK. C0MM0N 
BL0CK variables for the ADD code are given in Ref. 4. 

C0MM0N / CADD/ ADD code Air flow variables (Ref, 4) 

C0MM0N /CCRATE/ Variables for critical species rate equation 


ISDG 


Option flag 

SDG 

S 

8 ) 

Source of critical species 

AB.AF 

V A r 


EB, EF 

V E r 1 


ALPB.ALPF 

\ 

V° f 

1 Constants in Eq. 8.3.1 

BETB.BETF 

V 6 f ' 

i 

GAMB.GAMF 

"b^f 

) 

CARBN 

9 

c 

Carbon number 

R0C0NP 

R 

Gas constant (1.98717 cal/mole/°K) 

M0LE1.M0LE2 

V H 2 

Molecular weights 

Y02 

V 02 

Mole fraction of 0^ 

DF.LTB .DELTF 

{ b- 5 f 

Constants in Eq. 8.3.11 

ANUB.ANUF 

v b >v f 


ISGMAX 


Flag 

C0MM0N /CF/ 

Dependent Variable 


1-1 


F(J,K) 

f 

r v 

Dependent variable at 1-1, J,K 


1 

2 
I 
I 
I 
I 
I 

I 

a 

i 

rt 

II 
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C0MM0N / CFF/ 
A2(J) ,A3(J) 
B2(J) , B3(J) 
C(J) 

C0MM0N /CFG/ 
GF(J.K) 

C0MM0N /CFG2/ 
FG2(J,K) 


Coefficient of PDE 

A , A. Coefficient of second derivative 

2 3 

B^,B^ Coefficient of first derivative 

C Coefficient of function 

Dependent Variable 

f* Dependent variable at I,J,K 

J )N 

Mass fraction fuel 

C* Mass fraction of fuel 

J 


C0MM0N /CMU/ Turbulent viscosity 


XMUX(J) Turbulent viscosity (gm/cra/sec) 


C0MM0N /C0NV/ Conversion factors 


FTCM 

LBKG 

CKELRK 

FTM 


30.48 cm/ft 
.4538 lb/kg 
1.8 °R/°K 

.3048 m/ft 


C0MM0N /CT/ Coordinate transformation matrix 


TC0R(L,M,J,K) 


n L.M 


I.J.K 

C0MM0N /CX/ Cartesian coordinates 


Coordinate rotation matrix at I.J.K 


xa.J.K) 


Cartesian coordinate of point I,J,K 


C0MM0N /CY/ Computational coordinates at I 


Y(L.M.J.K) 

L - ] 

• 2 

*= 3 


Streamwise direction 
Principal normal direction 
Orthogonal direction 
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M - 1 

Y 

Coordinate at point I,J,K 

• 2 

h 

Metric (cm) 

- 3 

y 

Arc length (cm) 

C0MM0N /CYI/ 

Computational coordinates at 1-1 

YI(L,M, J,K) 


See C0MM0N /CY/ 

C0MM0N /DATABL/ I/O index 

data (see Table 9.1) 

INDX (N0P , 1) 


Unit number 

INDX (N0P,2) 


Block length 

INDX (N0P , 3) 


Pointer location 

INDX (N0P,4) 


Not used 

INDX (N0P,5) 



C0MM0N /RES1DL/ Residuals 


RES 


Residual 

IMAX 


Max no. iterations 

EPSL0N 


Tolerance 

PCT 


Percent change cutoff 

RLX 


Relaxation parameter 

C0MM0N / SAVC/ 

Solution variable on periodic boundary 

C2(J) 

C(J,1) 

Fuel mass fraction K«1 boundary 

CKLMl(J) 

C(J,KL) 

Fuel mass fraction K«1 boundary 

C0MM0N / SAVF/ 

Air flow variables 

FL0SAV(L, J) 


See ADD code 
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C0MM0N /S0URC/ 

Sources 

for fuel concentration equation 

W(J,K) 

V 

L 

3 

Source terns (gm/cm /sec) 

BDYY (L,M, J) 

• 

V 

LU 

2 

Wall source terms gm/cm / sec 

WPARM(L) 


See C0MM0N/MASSD/ in PTRAk 

C0MM0N/S0URC2/ 

Sources 

for critical species equation 


W2(J,K) 


Source for critical species equation 
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9. r > List of Local Variables 

A list of selected local variables appears below. These variables are 
alphabetically listed by subroutine name. 

ADDC0K (Ref. 4) 

i}(N,.l) ADD code variables 


1 LA 


No. streamwisc stations 

** 

it 

JLA 


No. normal stations 

•n 

I* 

1L.JL.KL 


No. mesh points i n Yj, Y>, Yj directions 

%} 

«# 

KADR 

r 

r 

Reference length 

mm 

ii 

DST1T 

As 

Y step size 

n 

Dl’lll 

•H 

Y step size 

i? 

AIN it; 



ti 




I! 

A1NT 


Area integral of f 

T] 

1 


Axial station 

: 1 

i f 

« 4 

I’CN 

l 

Integrand of integral 





cl 

KUNCT 


Durnrnv subrout Ine name 




• 

n 

:•{ 

C0KKK1 (Rel. ) 



• * 

r’( INVAR. 1 ,.n 


Air flow variables from ADD code 

• ♦ 

i 




.1 

SLUM 

Sc 

Schmidt number 

♦ t 

»| 

DSR 

11 

r 

Reference velocity 

« i 

R110K 

p 

r 

Re 1 e rence dens i t y 

- ? 

TKMl’R 

T 

r 

Reference t emp era lure 

r\ 

V 1 SLR 

u 

Re le rence viscosity 

t 

4 i 


i 
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PRFSR 

P 

r 

Reference pressure 

DSTEP 

As 

Y^ step size 

C0EFSG 



R0C0N 

R 

3 

Gas constant (82.0575 cm atm/mole/°K) 

R&C0NP 

R 

Gas constant (1.9817 cal/mole/°K) 

TINF 

T 

a 

Air temperature (°K) 

PJNF 

P 

a 

Air pressure (atm) 

UINF 

U 

Air velocity (m/sec) 

X02 

X Q2 

3 

Mole concentration of oxygen (mole/ cm ) 

Y02 

Y 02 

Mole fraction of oxygen 

DERBDY 



ID 


Direction index 

IV 


Variable index 

ITYPE 


Option 

DIN 

f 

Input array (dependent variable) 

11,12,13 


Dimensions of DIN 

YIN 

Y 

Input array (independent variable) 

JO,KO 


Smallest indices 

DG 

df/dy 

First derivative 

DDG 

2 2 
d f/dy 

Second derivative 

B,A 


Intermediate values 
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DERINT 

DERIV 

DFC0R 

WRATIO 

BIFFUS 

IADD 

IBEGIN 

IEND 

DM2DT 

DR0UT 

INUNIT 

ADDR 

BL0CK 

NM0VE 

FL0WIN 

C0NC C 

PBNDC 

AC,BC,CC,DC 

BDYVAL 

J ,K 

JJ,KK 

JBDY 

S0LN 


See subroutine DERBDY 
See subroutine DERBDY 

Fuel mass flow correction factor 

Add code strating station 
Starting station 
Ending station 
See C0MM0N/ CCRATE/ 

Unit number 
Output address 
Record length (words) 

Relative address 

Initial mass fraction of fuel 

Coefficient PDE 
Boundary value 
Mesh point 

Sub block mesh point 
Boundary index 
Solution 
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PERBC 

P0ISSN 

AC , BC , CC , DC 

JJ,XK 

S0L.N 

PRINPT 

PRRUN 

PRTDSK 

IDX1.IDX2 

ISTAT 

N0P 

RDINA 

I,J.K 

JJ.KK 

ADDR1 

N0P 

RDINB 

RESTRT 

I BEG IN 

PH0V 

RH 

RLXCAL 

RLXCAL 


See aubrouclne PBNDC 

Coefficients of PDE 
Sub block mesh point 
Solution 

See output format statements 
See C0MM0N/C0NV/ 

Indices of variables 

Uni t index 

Option 

Mesh point index 
Sub ^.ock index 
Block length 
See Table 9.1 
See subroutine RDINA 

Starting station no. 

2 

pu Mass flux (slug/ft /sec) / 

Relaxation parameter (Eq. 5.3.14) 
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S0LVE 

IADD 

IBEGIN , LEND 

J,K 

JJ.KK 

ISIDE 

FT1L 

RES 

S0LJRCE 

S 

SUMRY 

WR0UT 

WR0UTB 

WTFL0W 

WF 


f 


v+1 


w h 
L 1 


CPU 


ADD code station 
Beginning/Ending station 
Mesh point 

Sub block mesh point 
Boundary index 
Guess of solution 
Residual 

Integrand 

See output format statements 
See subroutine RDINA 
See subroutine RDINA 

2 

Mass flux fuel (slug/ft /sec) 
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